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Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

NONEQUILIBRIUM  FLUCTUATIONS  AND  TRANSPORT 
IN  SHEARED  FLUIDS 

BY 

JAMES  FRANCIS  LUTSKO 
December  1986 

Chairman:    James  W.  Dufty 
Major  Department:  Physics 

Hydrodynamic  fluctuations  in  a  sheared  fluid  far  from  equilibrium 
are  described  by  means  of  nonlinear  stochastic  differential 
equations.    Two  models  resulting  from  two  different  limits  of  the 
exact  Fokker-Planck  equation  describing  hydrodynamic  fluctuations  in 
an  arbitrary  fluid  are  presented.    The  first  model  results  from 
restricting  attention  to  fluctuations  which  take  place  on  length 
scales  large  with  respect  to  the  mean  free  path.    The  resulting 
Navier-Stokes-Langevin  equations  are  specialized  to  shear  flow  and  the 
associated  linear  hydrodynamic  modes  are  determined.    These  are  used 
to  calculate  the  correlation  functions  of  fluctuations  around  the 
steady  state  with  the  full  nonlinear  dependence  on  shear  rate 
necessary  for  large  shear  rates.    It  is  then  shown  how  the  coupling  of 
the  linear  modes  through  the  nonlinearities  in  the  stochastic 
equations  allows  one  to  calculate  the  effect  of  the  nonequilibrium 
state  on  the  macroscopic  transport  properties. 

vi 


The  second  model  results  from  the  short  time  limit  of  the  exact 
Fokker-Planck  equation  and  describes  fluctuations  even  on  atomic 
length  scales  provided  the  mean  free  path  is  small  enough.    The  linear 
modes  of  this  model  are  obtained  and  it  is  found  that  they  differ 
considerably  from  the  long  wavelength  modes.    In  particular,  it  is 
found  that  at  large  shear  there  are  two  sets  of  propagating  or  sound 
modes.    A  large  shear  instability  is  also  found  in  this  model  and  is 
tentatively  identified  with  a  large  shear  dynamic  phase  transition 
observed  in  computer  simulations  of  sheared  fluids. 
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CHAPTER  1 
INTRODUCTION 

Statistical  mechanics  is  a  methodology  for  determining  the 
physical  properties  of  many  body  systems  from  the  microscopic  dynamics 
of  its  constituents.    It  is  applicable  when  the  properties  one  wishes 
to  determine  arise  as  a  result  of  the  interaction  of  a  large  number  of 
the  elementary  constituents  or  when  the  time  scales  of  interest  are 
large  with  respect  to  the  typical  interaction  time  of  the 
constituents.    When  this  is  the  case,  attention  is  shifted  from  the 
dynamics  of  individual  constituents  to  the  density  in  phase  space 
(i.e.,  the  space  defined  by  the  degrees  of  freedom  of  the  system)  of 
the  system  as  a  whole  which  is  characterized  by  the  distribution 
function  of  the  system.    The  dynamics  of  the  distribution  function  is 
given  by  the  Liouville  equation,  which  expresses  the  conservation  of 
phase  space  density  in  a  conservative  system,  together  with 
appropriate  boundary  conditions,  initial  conditions  and  source  terms 
representing  the  action  of  nonconservative  external  interactions. 

When  the  distribution  function  is  constant  in  time,  the  system  is 
said  to  be  in  a  steady  state.    Equilibrium  is  the  particular  steady 
state  obtained  by  conservative  systems  and,  in  general,  is  unique  in 
that  the  associated  distribution  function  is  known.    For  this  reason, 
the  statistical  mechanics  of  equilibrium  systems  is  well  developed  and 
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quite  successful  at  predicting  macroscopic  properties  such  as 
transport  equations  and  thermodynamic  susceptibilities. 

The  statistical  mechanics  of  systems  not  in  equilibrium  is,  in 
contrast,  poorly  understood.    This  difficulty  arises  from  the  fact 
that  nonequilibrium  states  are  difficult  to  characterize;  the 
distribution  function  is  not  known.    For  systems  close  to  equilibrium 
linear  response  theory  provides  a  means  of  studying  nonequilibrium 
systems  by  relating  their  properties  to  those  of  equilibrium  systems 
(e.g.,  fluctuation-dissipation  theorems  relate  transport  properties  to 
equilibrium  fluctuations).    However,  far  from  equilibrium  one  has  only 
the  Liouville  equation. 

To  get  a  handle  on  systems  far  from  equilibrium,  one  typically 

restricts  attention  to  some  subset  of  the  degrees  of  freedom  which 

obeys  an  approximately  closed  dynamics  on  the  length  and  time  scales 

of  interest.    Closed,  in  this  context,,  means  that  couplings  to  the 

degrees  of  freedom  not  in  the  subset  being  studied  are  negligible  on 

the  average.    For  example,  the  macroscopic  equations  of  continuum 

mechanics  for  the  local  conserved  densities  (e.g.,  hydrodynamics) 

often  provide  an  adequate  characterization  of  the  nonequilibrium 

state.    The  problem  then  becomes  one  of  determining  the  average 

dynamics  of  these  variables  (i.e.,  macroscopic  transport  equations) 

and  their  fluctuations.    At  present,  almost  nothing  is  known  about  the 

transport  properties  of  systems  far  from  equilibrium. 

This  semi-phenomenological  approach  was  first  proposed  some  fifty 

years  ago  by  Onsager1  and  has  essentially  been  verified  for  systems 
2  3 

near  equilibrium.  '      Recently,  it  has  been  extended  to  apply  to 
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systems  far  from  equilibrium,  ' D  but  it  is  difficult  to  evaluate  the 
success  of  such  an  extension  since  it  has  been  applied  to  very  few 
problems.    The  reason  it  has  not  been  more  thoroughly  tested  is  that 
its  application  is  straightforward  only  if  the  nonequil ibrium  state  is 
characterized  by  a  few  simple  macroscopic  properties. 

One  such  nonequilibrium  system  is  a  simple  fluid  in  a  state  of 
uniform  shear  flow.    In  uniform  shear  flow,  the  pressure  and 
temperature  are  spatially  constant  while  the  flow  velocity  varies 
linearly  with  distance  in  a  direction  perpendicular  to  the  flow 
velocity  (Figure  1.1).    In  a  rectangular  coordinate  system  (denoted  as 
x-y-z),  the  flow  field,  u^r.) ,  may  be  written 

u^tr)  =  a-r  (1.1) 

where  the  shear  tensor,  a,  is 

a  =  axy  (1.2) 

and  "a"  is  the  shear  rate.    The  deviation  of  this  system  from 
equilibrium  is  controlled  by  a  single  parameter,  the  shear  rate, 
giving  a  simple  characterization  of  the  nonequilibrium  state. 

A  sheared  fluid  is  one  of  the  simplest  nontrivial  nonequilibrium 
systems  imaginable  and,  indeed,  may  appear  so  simple  as  to  be 
uninteresting.    However,  even  in  this  simple  system  drastic  effects 
due  to  deviations  from  equilibrium  are  present,  for  example  in  its 
transport  properties,  giving  the  system  an  intrinsic  as  well  as 
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pedagogical  interest.    In  the  remainder  of  this  introduction,  we  will 
first  describe  three  anomalies  in  the  behavior  of  sheared  fluids,  and 
then  we  will  describe  the  methods  and  physical  ideas  which  will  be 
used  to  attempt  an  understanding  of  the  system.    The  introduction 
concludes  with  an  outline  of  the  work  to  be  presented. 

Sheared  fluids  are  of  value  as  model  nonequilibrium  systems  not 
only  because  of  the  simplicity  of  the  macroscopic  state,  but  also 
because  a  great  deal  of  effort  has  gone  into  developing  methods  of 
simulating  such  systems  on  a  computer . 6-1 ^    This  is  of  prime 
importance  because,  as  will  be  discussed  in  detail  later,  laboratory 
experiments  involving  simple  atomic  fluids  are  at  present  limited  to 
the  small  shear  rate  regime  (i.e.,  close  to  equilibrium).    As  will  be 
discussed  in  the  next  chapter,  experiments  are  possible  involving  more 
complex  rheological  fluids,  such  as  polymeric  fluids  and  glasses,14 
and  colloidal  suspensions1 5  which  exhibit  many  of  the  properties  to  be 
discussed  here  for  simple  fluids.    We  concentrate  on  simple  fluids, 
however,  to  simplify  the  theoretical  description.    Thus,  the  only 
"empirical"  data  which  exist  for  simple  fluids  at  large  shear  rates 
are  from  simulations.    In  fact,  it  is  this  body  of  "empirical"  results 
which  has  spurred  most  of  the  recent  theoretical  interest.1*'16"25 

The  first  important  anomaly  observed  in  these  simulations  was 
associated  with  the  phenomena  of  shear  thinning.    As  some  fluids  are 
sheared,  their  viscosity  drops;  they  become  "thinner"  or  less 
viscous.    The  dependence  of  viscosity  on  shear  rate  had  been 
calculated  in  1969  by  Yamada  and  Kawasaki16  as 
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n(a)  =  nQ  -  ^a172  +  0(a).  (1 .3) 

The  nonanalytic  nature  of  the  dependence  on  the  shear  rate  was 
previously  unsuspected  and  indicates  the  breakdown  of  nonlinear 
response  based  on  a  formal  expansion  of  the  Liouville  equation  in  the 
shear  rate.    In  addition,  the  result  of  the  simulations  was  that  ni 
was  measured  to  be  two  orders  of  magnitude  larger  than  Yamada  and 

16  1  ft 

Kawasaki      predicted.    Several  other  authors      calculated  rii  and 

obtained  different  numbers,  although  all  are  two  orders  of  magnitude 

too  small.    A  recent  proposal  by  van  Beijeren      for  the  resolution  of 

this  problem  will  be  discussed  below. 

In  1982,  a  second  anomaly  observed  in  simulations  was  reported  by 
1  2 

Erpenbeck.       He  found  that,  at  very  high  shear  rates,  the  system 
passed  from  the  amorphous  fluid  phase  into  an  ordered  phase  in  which 
the  particles  line  up  one  behind  another  in  hexagonally  packed 
"strings"  (Figure  1.2).    There  has  recently27  been  some  debate  about 
the  actual  nature  of  the  string  phase  but  it  is  generally  agreed  that 
a  high  shear  rate  instability,  or  dynamical  phase  transition, 
exists.    We  shall  present  below  a  model  exhibiting  an  instability  on 
atomic  scales  which  matches  up  qualitatively  quite  well  with  that  seen 
in  the  simulations  (note  that  macroscopic  hydrodynamics  is  believed  to 
be  stable  for  planar  shear  flow28). 

A  third  anomaly  has  been  reported  by  Evans7  as  occurring  in  two- 
dimensional  fluids.    His  simulations  show  the  shear  viscosity 
remaining  constant  as  the  shear  rate  increases  until  at  a  critical 
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shear  rate  the  fluid  begins  to  thin  in  the  expected  manner.  Evans 
attributes  this  to  a  small  shear  rate  instability  for  which  no 
theoretical  understanding  exists  as  yet.    However,  the  transport 
properties  of  two  dimensional  fluids  are  anomolous  even  in  equilibrium 
(two  is  the  critical  dimension)  leading  one  to  believe  that  Evans' 
observations  may  be  related  to  this  fact. 

Most  of  the  work  to  be  presented  here  will  consist  of  an  attempt 
to  construct  and  analyze  models  of  fluctuations  around  the 
nonequilibrium  state  in  dense  sheared  fluids.    Once  the  fluctuations 
are  statistically  characterized,  one  is  able  to  calculate  time 
dependent  correlation  functions  which  in  turn  can  be  related  to 
observable  physical  properties  of  the  system.    The  phenomenological 
ideas  behind  this  approach  are  very  similar  to  those  in  the  Einstein- 
Smolochowski  model  of  Brownian  motion. 2^    One  considers  the  dynamics 
of  the  hydrodynamic  variables  (analogous  to  the  position  of  the 
Brownian  particle)  around  the  steady  state.    The  hydrodynamic 
variables  are  appropriate  to  describe  properties  on  a  scale  which  is 
large  compared  to  the  microscopic  correlation  length.    Because  these 
variables  are  the  local  conserved  densities,  it  is  assumed  that  they 
decay  according  to  the  macroscopic  conservation  laws  (e.g.,  their  time 
evolution  is  governed  by  the  Navier-Stokes  equations)  with  the 
exception  that  the  equations  include  stochastic  sources  which 
represent  the  thermal  noise  responsible  for  the  creation  of  the 
fluctuations.4'30    Finally,  the  stochastic  sources  are  generally  taken 
to  be  Gaussian  white  noise  with  correlations  fixed  by  a  fluctuation- 
dissipation  theorem. 
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Having  these  equations,  one  in  principle  solves  them  for  the 
hydrodynamic  variables  as  functionals  of  the  stochastic  sources  and 
thus  may  calculate  time  dependent  correlation  functions  using  the 
known  statistical  properties  of  the  sources.    As  alluded  to  above, 
these  correlation  functions  can  in  turn  be  used  to  calculate  various 
physical  quantities  such  as  the  effect  of  the  steady  state  on 
transport  coefficients  and  thermodynamic  susceptibilities. 

While  providing  a  closed  and  intuitively  appealing  model  of 
fluctuations,  this  description  is  entirely  phenomenological  and  as 
such  represents  an  uncontrolled  approximation.    Several  authors5'^1 
have  addressed  the  question  of  the  formal  derivation  and  validity  of 
fluctuating  hydrodynamics  and  we  shall  outline  their  ideas  in  a  later 
chapter.    The  result  of  these  studies  is  that  the  foundations  and 
validity  of  this  approach  are  much  better  understood. 

Thus,  one  approach  to  the  study  of  nonequilibrium  systems  is  to 
replace  the  microscopic  description  of  the  system,  as  characterized  by 
the  distribution  function  and  Liouville's  equation,  with  a 
hydrodynamic  description.    A  hydrodynamic  description  of  a  system  is 
obtained  as  a  result  of  two  approximations.    The  first  is  that 
attention  is  restricted  to  the  local  densities  of  the  conserved 
variables  (e.g.,  mass,  momentum  and  energy  in  a  simple  fluid).  On 
length  scales  which  are  large  compared  to  the  microscopic  correlation 
length  (i.e.,  the  mean  free  path,  imfp) ,  these  variables  are  expected 
to  decay  slowly  since  their  integrals  over  the  whole  system  do  not 
change  at  all  with  time.    Thus,  on  these  length  scales,  these 
variables  are  assumed  to  approximately  decouple  from  the  remaining, 
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quickly  fluctuating,  microscopic  degrees  of  freedom  which  give  rise  to 
dissipation  and  thermal  fluctuation  of  the  slow  variables.  Concurrent 
with  this,  if  k  is  the  wavevector  (inverse  wavelength  or  gradient  in 
real  space)  of  the  conserved  densities,  it  is  required  that  (k£mfp)  be 
small.    The  equations  governing  the  dynamics  of  the  conserved 
densities  may  then  be  expanded  in  this  small  parameter  with  truncation 
at  second  order  yielding,  for  a  fluid,  the  usual  Navier-Stokes 
equations . 

In  a  sufficiently  dense  fluid,  the  mean  free  path  is  actually 
smaller  than  the  typical  size,  a,  of  the  atoms  making  up  the  fluid. 
In  this  case,  the  hydrodynamic  description  will  include  terms  of  all 
order  in  (ka)  while  still  expanding  in  (k&mfp).    The  result  is  a 
hydrodynamic  description  of  the  (dense)  fluid  which  contains 
information  about  the  structure  of  the  fluid  and  is  valid  on  atomic 
(i.e.,  ko  -  1)  length  scales  although  still  restricted  to  lengths 
which  are  large  relative  to  the  mean  free  path  (ki.mfp  <<  1).  One 
might  think,  and  indeed  it  was  long  thought,  that  atomic  scale 
fluctuations  require  knowledge  of  all  small  scale  fluctuations;  after 
all,  the  reason  for  considering  only  conserved  densities  at  long 
wavelengths  is  that  they  decay  on  macroscopic  length  and  time  scales 
while  all  other  fluctuations  decay  on  much  shorter,  microscopic, 
scales.    However,  while  studying  models  of  small  length  scale 
fluctuations,  de  Schepper  and  Cohen-5    discovered  that  even  in  this 
regime  the  conserved  densities  decayed  much  slower  than  other 
variables  and  gave  the  dominant  contribution,  on  these  scales,  to  the 
density-density  correlation  function.    This  indicated  that  restricting 
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attention  to  the  hydrodynamic  subspace  was  a  valid  means  of  studying 
fluctuations  even  on  very  small  length  scales.    Thus,  hydrodynamic 
models  can  be  used  to  study  even  atomic  length  scale  fluctuations 
provided  the  mean  free  path  is  small  enough.    In  particular, 
van  Beijeren's  explanation     of  the  discrepancy  between  theory  and 
simulation  results  for  the  size  of  the  shear  rate  dependent 
renormalization  of  the  shear  viscosity  is  based  on  this  generalization 
of  fluctuating  hydrodynamics  to  include  atomic  scale  fluctuations  in 
dense  fluids. 

The  work  to  be  presented  below  consists  of  modeling  hydrodynamic 
fluctuations  as  a  means  of  understanding  the  nonequilibrium  state.  In 
Chapter  2,  we  will  discuss  experimental  and  simulational  means  of 
studying  steady  state  shear  flow  and  the  connection  between' the  two. 
The  third  chapter  discusses  the  formal  ideas  of  statistical  mechanics 
used  to  construct  models  of  nonequilibrium  fluctuations  and,  in 
particular,  to  construct  fluctuating  hydrodynamics  at  long  wavelengths 
as  well  as  generalized  fluctuating  hydrodynamics.    This  is  followed  by 
Chapters  4  and  5  analyzing  the  long  wavelength  model,  the  linearized 
and  nonlinear  versions  respectively,  with  the  aim  of  deriving  various 
physically  observable  consequences  of  the  nonequilibrium  state.  The 
sixth  and  seventh  chapters  contain,  respectively,  the  development  and 
analysis  of  the  small  wavelength  model  and  include  a  discussion  of  the 
hydrodynamics  and  large  shear  instability  occurring  on  these  length 
scales.    The  dissertation  concludes  with  a  discussion  of  the  results 
obtained,  the  completeness  of  the  present  description  and  future 
directions . 


CHAPTER  2 
PHENOMENOLOGY  AND  SIMULATION 

In  this  chapter,  we  discuss  the  basic  phenomenology,  such  as 
length  and  time  scales  and  strength  of  shear,  which  affects  the 
observability  of  deviations  from  equilibrium  behavior.    We  then 
discuss  the  possibility  of  performing  laboratory  experiments  which 
measure  nonequilibrium  properties.    Finally,  we  discuss  the  methods  of 
simulating  uniform  shear  flow  on  a  computer. 

Phenomenology 

In  macroscopic  hydrodynamics,  one  typically  has  three  types  of 

parameters:    thermodynamic  derivatives  such  as  the  sound  velocity,  c0, 

transport  coefficients  such  as  the  kinematic  viscosity,  v0,  and  the 

wavevector,  k,  the  magnitude  of  which  is  inversely  proportional  to  the 

length  scale  of  the  phenomena  one  is  interested  in.    The  time  scale  of 

2  -1 

hydrodynamic  dissipation  £s  (vQk  )      (in  equilibrium,  momentum  fluctu- 

~v0k  t 

ations  dissipate  as  e  ) .    In  hydrodynamics,  one  usually  uses 

perturbation  theory  treating  the  wavevector  as  a  small  parameter.  The 
dimensionless  parameter  which  actually  occurs  is  (v0k/cQ)  thus 
requiring  that  vQk  <<  cQ  or  vQk    <<  cQk.    Physically,  this  means  that 
the  rate  of  hydrodynamic  dissipation  is  small  compared  to  the  rate  of 
sound  propagation.    This  is  consistent  with  the  Navier-Stoke' s 
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equations  governing  hydrodynamics  since  they  are  themselves  the  result 

of  an  expansion  to  second  order  in  the  parameter  (ki,      ),  where 

mfp 

amfn  is  the  mean  free  Path«  and  in  fact  cn/vn  ""  *  -  •    Typical  values 
miP  0    0       mfp  2 

of  these  constants  for  water  at  STP35  are  v    -  1  .8  •  10~2 

0  s 

5  -7 
0Q  -  1.5  •  10    cm/s  so  I       =  10      cm.    Thus,  the  applicability  of 

hydrodynamics  requires  that  k  <<  8.3  •  1 06  cm"1 ,  or  ko  <<  10~1 
~8 

where  a  -  10      cm  is  a  typical  atomic  radius. 

When  a  fluid  is  sheared,  a  new  time  scale,  the  inverse  of  the 

shear  rate  (denoted  here  as  "a"),  is  introduced.    When  the  shear  rate 

is  as  large  as  the  rate  of  dissipation,  we  expect  that  large 

deviations  from  equilibrium  behavior  will  occur.    For  this  reason,  we 
2 

refer  to  a  -  vQk    as  the  large  shear  regime.    When  studying  the  long 

p 

wavelength  hydrodynamic  model,  we  will  require  that  a  <  vQk  thus 
including  the  large  shear  regime  while  allowing  us  to  use  perturbation 
theory  (since,  given  this  bound,  a  <<  cQk). 

Experiments 

Various  methods  exist  for  measuring  the  static  and  transport 
properties  of  fluids.    For  example,  neutron  scattering  can  be  used  to 
measure  the  static  and  dynamic  structure  factors  (related  to  the  pair 
correlation  function  and  transport  properties,  respectively)  on  atomic 
length  scales.       Light  scattering  samples  the  dynamic  structure 
factor  on  hydrodynamic  length  and  time  scales.37    These  methods  and 
others  are  available  to  measure  the  corresponding  properties  for 
nonequilibrium  fluids.    However,  the  present  ability  to  test  theories 
of  sheared  systems  is  limited  by  the  shear  rates  which  can  be  obtained 


in  the  laboratory.    At  present,  the  largest  shear  rates  reported  are 

less  than  10    HZ.^0    This  shear  rate  is  in  the  large  shear  regime 

1/2  2-1 
when  k  -  (a/vQ)       -  7.45  •  10    cm    .    This  means  that  until  we  look 

at  length  scales  as  large  as  about  10     cm  the  effects  of  shear  will 

be  small  compared  to  the  usual  equilibrium  properties.    Thus,  light 

scattering,  which  measures  properties  on  a  length  scale  of  the 

-5 

incident  beam's  wavelength  -  10      cm,  and  neutron  scattering,  which  is 

Q 

not  practical  for  wavelengths  above  about  10  °  cm,  require  shear  rates 
about  two  orders  of  magnitude  larger  than  those  attainable  in  the 
laboratory. 

There  are  systems,  other  than  simple  fluids,  in  which  similar 

nonequilibrium  effects  occur  and  which  have  the  practical  advantage 

that  they  have  long  relaxation  times.    Because  the  large  shear  regime 

p 

occurs  for  a  -  inverse  relaxation  time  (-  vQk  in  a  simple  fluid)  it 

thus  becomes  accessible  in  these  systems.    One  example  of  such  a 

system  is  a  colloidal  suspension  in  which  the  relaxation 
2 

time  t  -  d  /D  for  a  typical  interparticle  spacing,  d,  and  diffusion 
constant  of  the  spheres,  D.    In  this  system,  the  density  of  spheres 
determines  the  relaxation  time  and  the  nature  of  the  probe  is  not 
restricted  as  in  a  simple  fluid.    Experimental  studies  of  sheared 
colloidal  suspensions  have  been  performed1^  and  numerous 
nonequilibrium  phenomena  are  observed  including  a  string  phase  similar 
to  that  seen  in  simulations  of  simple  fluids. 

Another  class  of  systems  exhibiting  easily  accessible 
nonequilibrium  phenomena  is  macromolecular  (e.g.,  polymeric)  fluids 
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and  glasses.       The  polymeric  systems  are  not  simple  fluids  because 
they  are  composed  of  anisotropic  molecules.    However,  Simons  et  al.^ 
have  shown  that  strong  similarities  exist  between  the  phenomena  seen 
in  simple  and  complex  fluids.    These  fluids  demonstrate  a  wide  variety 
of  non-Newtonian  behaviors  including  shear  thinning.  The 
accessibility  of  these  phenomena  again  arise  due  to  the  fact  that  the 
constituents,  the  macromolecules,  have  long  relaxation  times  lowering 
the  lower  bound  of  the  large  shear  regime. 

In  this  work,  we  will  restrict  our  attention  to  simple  fluids. 
As  the  preceding  discussion  indicated,  we  will  not  have  direct  experi- 
mental tests  of  our  results.    The  reason  we  do  not  consider  colloidal 
suspensions,  polymeric  fluids  or  some  other  such  system  is  that  the 
results  from  computer  simulations  of  simple  fluids  have  shown  that 
even  our  understanding  of  the  simplest  system  is  incomplete.  However, 
because  ef fects,  such  as  shear  thinning  and  the  string  phase,  seen  in 
the  simulations  are  observed  in  physically  related  laboratory  systems, 
we  have  some  confidence  in  the  simulations.    Thus,  we  propose  to 
attempt  to  understand  the  simplest  sheared  system  and  hope,  in  the 
process,  to  gain  some  understanding  by  analogy  of  the  more  complicated 
systems  in  which  the  nonequilibrium  state  has  important  physical 
consequences. 

Simulations 

Most  of  the  data  for  simple  fluids  have  been  generated  by 
nonequilibrium  molecular  dynamics6-1 3  rather  than  in  the  laboratory. 
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Such  simulations  require  both  the  establishment  of  the  linear  shear 
profile  as  well  as  the  maintenance  of  constant  temperature. 

In  equilibrium  molecular  dynamics,  one  considers  a  box  with 
periodic  boundary  conditions  containing  N  particles.    At  each  time 
step  in  the  simulation,  one  solves  the  equations  of  motion  (e.g., 
Newton's  equations)  for  all  of  the  particles  and  moves  them 
accordingly.    The  periodic  boundary  conditions  have  the  effect  of 
decreasing  the  size  (number)  dependence  of  the  results  allowing  one  to 
achieve  the  bulk,  or  thermodynamic,  limit  with  on  the  order  of  10^ 
particles.    In  this  way,  one  may  study  static  properties,  such  as  the 
pair  correlation  function,  and  dynamic  properties  like  the  decay  of 
correlations.    Clearly,  simulations  have  the  advantage  over 
experiments  in  that  one  can  obtain  much  more  data  about  the  system 
since  one  has  the  microscopic  state  at  one's  disposal.    They  are, 
however,  limited  by  the  number  of  particles,  N,  and  the  length  of  time 
which  can  be  simulated.    In  particular,  the  size  restriction  limits 
the  smallest  wavevectors  which  can  be  studied. 

The  linear  shear  profile  is  generated  by  modifying  the  simple 
periodic  boundary  conditions  of  equilibrium  simulations.    The  system 
remains  periodic  in  the  x  and  z  directions  but  when  a  particle  with 
momentum  £  passes  through  the  y  =  L  plane  at  time  t  and  with  coordi- 
nates (x,L,z),  it  is  reentered  at  y  =  0  with  coordinates  (x-aLt,  y,  z) 
and  momentum  £  -  m§Lx.    These  boundary  conditions  are  known  as  Lee- 
Edwards'  boundary  conditions.    They  become  more  transparent  when 
expressed  in  the  local  rest  frame  coordinates,  q  and  jj    defined  as 
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q  =  q  -  a  •  £t 


£  -  ma  •  £. 


(2.1) 


In  these  coordinates,  the  local  flow  velocity  is  zero  and  it  is  seen 
that  the  Lee-Edwards'  conditions  correspond  to  simple  periodic 
boundary  conditions  in  the  rest  frame. 

The  temperature,  T,  is  a  measure  of  energy  fluctuations  around 
the  steady  state  and  is  defined  by 

3  1    N  2 

2  NkBT "  <^  l  (£i  "  mi  '  aj)  > 

where  kg  is  Boltzmann's  constant.    Substitution  of  the  conditions  of 
spatially  uniform  shear  flow  into  the  Navier-Stokes  equations  shows 
that  a  sheared  fluid  will  undergo  viscous  heating.    Specifically,  one 
finds  that 

3e  2 

^  -  na  ,  (2.3) 
where  e  is  the  internal  energy  density,  leading  to 


c(t)  =  e(0)  +  na2t. 


(2.4) 
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For  water  sheared  at  10    Hz  at  STP  this  leads  to  a  heating  rate  of 
0.024  K/s.    For  temporally  uniform  shear  flow,  particularly  for  the 
very  high  shears  applied  in  computer  simulations,  an  external, 
nonconservative  force  is  necessary  to  hold  the  temperature  constant. 
This  heat  extracting  interaction  is  called  a  thermostat.    This  is 
equivalent  to  adding  a  drag  force  to  the  equations  of  motion: 

d 

HtRi        "  A(r)  £i 

A(D  =  I  L  •  (F    -  a  •  £  )  /  £  g2.  (2.5) 
J  J  J        j  J 

where  the  form  of  A(D  has  been  chosen  to  hold  the  temperature  as 
defined  in  eq.  (2.2)  constant. 

The  equation  of  motion  form  for  the  thermostat  is  important 
because,  with  the  Lee-Edwards*  boundary  conditions  for  imposing  the 
flow  field,  the  simulations  can  be  mapped  onto  a  well  defined 
mathematical  problem:    the  solution  of  Liouville's  equation,  modified 
to  account  for  the  thermostat,  with  the  imposed  boundary  conditions. 
The  modified  Liouville  equation  is  obtained  from  the  general  form  of 
Liouville's  theorem,^ 

c!t"  +  iUi  '  ai~  +  h  '  ff^ptr.t)  -  aj:  •  CA(r)£.p(r;t)]  =  o 

(2.6) 
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where  the  repeated  index,  i,  is  to  be  summed  and  runs  from  1  through  N 
and  the  dot  indicates  an  ordinary  vector  scalar  product.    So,  (2.6) 
may  be  rewritten  as 


ft  P  +  £p  =  0,  (2.7) 


with 


£  =  £n  +  £.  .   +  £    .  (2.8) 
0       int  ext 


where  the  free  streaming  and  interaction  terms  are,  respectively, 


£i„t  -  h  ■  4  <2-9> 


while  the  external  term,  due  to  the  nonconservati ve  thermostat,  acts 
on  the  distribution  function  as 


*ext  p(F;t)  =  "  9i~  '  ^(Dfi^r.t)].  (2.10) 


Because  this  operator  includes  the  nonconservation  of  volume  in  phase 

space  due  to  the  nonconservative  force  (Liouville's  theorem  states 

tnat       (p(r)dr)  =  0  and  with  nonconservative  forces  —  dr  *  0) 
Qt  dt 

It  has  a  different  form  when  acting  on  any  other  phase  function. 
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Specifically,  if  A(r)  is  an  arbitrary  phase  function,  we  have  that 


The  bar  indicates  the  difference  between  these  operators  and  we  note 
that  £Q  =  £Q.    It  is  easy  to  see  that  £  and  £  are  adjoints  with 
respect  to  phase  space  integration. 

Finally,  we  should  mention  that  when  working  with  hard  spheres, 
the  potential  is  singular  and  the  interparticle  force,  appearing 
in  these  equations  ill  defined.    This  difficulty  is  eliminated  by 
replacing  the  interaction  term  with  the  so-called  T-operators  which 
generate  the  singular  hard  sphere  dynamics.    Henceforth,  for  hard 
spheres,  we  will  write 


where  ±  refers  to  forward  and  backward  time  equations.  The  form  of 
these  operators  can  be  found  in  the  literature110  and,  for  example, 


£extA(D  -  -  Kr)  £.  •  Jj-A(r). 


i 


(2.11) 


V^i'W 


(2.12) 


a  • 


m 


6<a12  -  r) 


•  Cb12  -  1]  A(xlfx2)} 


(2.13) 


where  x1  =      .p^  ,  etc.;  a  is  the  hard  sphere  diameter,  p_. 


•12  =  £1  -  £2 


and 


b12A(£i  •£!  ;a2.£2)  =  A^.Et  "       •  £1;a2,£2  +  oo  •  £12).  (2. 


We  note  here  that  the  T-operators  are  not  self  adjoint,  as  is  the 
interaction  term  for  continuous  potentials,  eq.  (2.9),  so  their 
adjoints  are  denoted  with  bars  as  well. 


CHAPTER  3 
FLUCTUATION  MODELS 

In  this  chapter,  we  discuss  the  statistical  mechanical  basis  of 
the  models  to  be  analyzed.    Specifically,  the  models  we  will  study 
will  be  seen  to  be  two  different  limits  of  a  general  Fokker-Planck 
equation.    The  detailed  form  of  the  models  will  also  be  given. 

Generalized  Langevin  Equation 
Our  basic  approach,  as  discussed  in  the  introduction,  to  studying 
nonequilibrium  fluctuations  around  shear  flow  is  to  choose  a 
restricted  set  of  variables  whose  dynamics  we  wish  to  study.    At  a 
fundamental  level,  any  classical  system  of  N  structureless  particles, 
or  atoms,  is  described  by  the  6N  coordinates,  collectively  called  r, 
which  are  the  positions,  q[t  and  momenta  p_    of  the  atoms  (where  i 
ranges  from  1  to  N).    The  distribution  function  of  the  system, 
p(T;t),  which  describes  the  density  of  states  in  phase  space, 
satisfies  the  Liouville  equation, 

|jT  p(r;t)  +  £  P(r;t>  =  0;  (3.1) 

the  explicit  form  of  which  was  discussed  in  a  Chapter  2.    The  set  of 
restricted  variables  we  choose  to  study  is  the  set  of  conserved 
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densities.    These  are  mass  density,  p(r),  energy  density,  e(r), 

and  momentum  density,  j  (r).    A  caret  over  a  function  indicates  that 

it  is  a  phase  function.    Their  definitions  are 

N 

p(r)  =    I    m  A(r-q  ) 
~       n-1         ~  ~n 

N  N 
e(£)  =    *         (rJ    A(r-a  )  +  ±    I      V(  q  -q  ,|)  A(r-q  ) 
n=1    m     n  n       2  n*n*         n    n  n 

N 

Ja(n)  -    Z    £n  Mr-a  )  (3.2) 
n=1  n 

where  A(r)  is  some  localized  function  of  r  (for  instance,  a  Dirac 
delta  function).    These  variables  will  be  denoted  collectively  as 
y(r).    Our  reasons  for  choosing  this  particular  set  were  alluded  to  in 
the  introduction.    At  length  scales  large  with  respect  to  the  mean 
free  path,  we  expect  that  the  conserved  densities  are  the  only 
variables  with  long  relaxation  times  (i.e.,  long  with  respect  to  the 
collison  frequency).    Therefore,  we  expect  that  the  statistical 
properties  of  the  system  at  these  scales  can  be  adequately  described 
in  terms  of  these  variables  with  the  other  degrees  of  freedom  being 
modelled  as  rapidly  varying  sources.    In  dense  fluids,  the  mean  free 
path  can  be  smaller  than  the  atomic  dimensions  and  it  is  found,  as  one 
might  guess,  that  the  hydrodynamic  variables,  or  conserved  densities, 
still  have  long  relaxation  times  and  so  still,  on  the  average,  obey  an 
approximately  closed  dynamics.    Because  the  other  degrees  of  freedom 
are  more  rapidly  varying,  they  may  again  be  modelled  as  stochastic 
sources. 
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In  the  analysis  to  follow,  it  will  prove  convenient  to  work  with 
the  Fourier  transform  of  the  variables  in  eq.  (3.2).  This  is  defined 
to  be 


ik  t 

y(k)  =  V"1  /  dr  e  "  ~  y(r).  (3-3) 


We  also  note  that  the  vector  y(k)  is  explicitly, 


y(k)  =  (p(k),  i(k),  ji(k) ) .  (3.4) 


For  notational  convenience,  we  will  write 


ya  **  yct(-)  (3*5) 


where  the  Latin  index  includes  both  the  vector  index  a  and  the 
wavevectorJ<_.    Repeated  Latin  indices  are  to  be  summed.    The  average 
of  ya.  or  ya(k),  over  the  nonequilibrium  ensemble,  p(T;t),  will  be 


denoted  as  <y  ;t>  where 

Si 


<ya;t>  =  /  dr  P(r;t)ya  (3.6) 

and  are  the  macroscopic  variables  one  would  measure  in  the  laboratory. 

In  general,  one  also  defines  time  dependent  phase  functions, 
y(k;t),  as 
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where 


y  (k;t)  =  e"£t  y  (k) 
a  —  a  — 


ya(t)  (3.7) 


<ya;t>  =  /dr  ya  e"£t  P(r;o) 


=  /dr  [e'£t  y  ]  p(r;0) 

CL 


=  /dr  y  (t)  P(r;o) 


<y.(t);0>.  (3.8) 

a 


So,  the  macroscopic  variables  are  the  average  over  an  initial 
ensemble,  p(T;0),  of  y  (t).    The  fluctuations  of  y  (t), 
denoted  as  6y  (t),  are  defined  to  be 

a 


6ya(t)  =  ya(t)  "  <Vt)!0>-  (3.9) 


Because  we  are  restricting  attention  to  the  set  of  variables 
y(t),  it  is  useful  to  introduce  the  distribution  function  of  these 
variables,  p(y;t),  defined  to  be 
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p(y;t)  =  Jar  5(y-y')  p(r';t),  (3.10) 

where  y  is  some  possible  value  of  the  phase  function  y.    So,  for  any 
function  of  the  y's,  say  A[y],  it  is  clear  that 

<A[y];t>  =  Jdr  A[y]  p(T;t) 

-  Jdy  A[y]  p(y;t).  (3.11) 


Finally,  for  a  steady  state  such  as  uniform  shear  flow,  we  have 
that  Pag(r;t)  =  Pgs(r;0)  =  P    (T).  So, 


<y  ;t>      =  <y  ;0> 
a     ss       •'a  ss 


sy0a.  (3.12) 


Grabert  et  al.  1  show,  using  the  projection  operator  technique 

hp 

introduced  by  Zwanzig,      that  the  distribution  p(y;t) 

satisfies  the  following  exact  generalized  Fokker-Planck  equation, 


3  9  t 

3t  P(y;t)  +  —  {vjy]  p(y;t)  -  JQds  Jdy'  D^Cy.y'js]  Pss(y') 


3y. 


3yb'  pss(y') 
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where  p    (y)  is  the  steady  state  distribution,  v  [y]  is  called  the 

S3  £1 

drift  vector  and  D^Cy.y'js]  is  called  the  diffusion  matrix.  Repeated 
Latin  indices  are  again  to  be  summed  in  (3.13).    The  functionals  v  and 
D  are  expressed  in  terms  of  steady  state  averages  and  their  detailed 
form  is  given  in  Grabert  et  al.  1    Here,  we  only  note  that 


—  {v  [y]  P  (y)}  =  0  (3.14) 
3ya 

so  P    (y)  explicitly  satisfies  (3.13). 

3  3 

Associated  with  (3.13)  is  an  exact  generalized  Langevin  equation 
for  the  variables  y  (t), 


3t 


ya(t)  -  vaCy(t)]  +  /Jds  {Dab[y(t-s);s]  Fb[y(t-s)] 


3yb(t-s) 


Dab[y(t-s);s]} 


Ra(t)  (3.15) 


where  the  thermodynamic  forces,  Fb[y(t-s)]  are  defined  in  terms  of  the 
steady  state  distribution  as 

3  in  P  (y) 

V'*-"  ly-J(t-s)  o.«) 

which  also  appears  in  (3.13).    The  functional  va  is  the  same  as  that 
appearing  in  (3.13)  while 
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Dab[y;s]  =  /dy'  DabCy,y  ;s]|y_^  (3.17) 


The  sources  R  (t)  have  the  following  constrained  averages, 

a 


0  =  /dl"  6(ya-ya')  Ra'(t)  Pga(r») 


Dab[y:t]  =  /dI"  6(Vya')  Ra'(t)  Rb'(0)  pss(r,)  (3'18) 

where  the  hat  on  Ra(t)  indicates  it  to  be  a  phase  function  and  where  a 

primed  function  has  argument  I".    The  first  of  equations  (3.18) 

indicates  that  Ra,(t)  is  orthogonal,  under  the  steady  state  average, 

to  any  function  of  y  and  indicates  that  the  constrained  average  of 

R  (t)  at  fixed  y  is  zero.    The  second  of  equations  (3.18)  has  the  form 

of  a  fluctuation-dissipation  theorem  relating  the  correlation  of  the 

"fluctuating  force"  to  the  diffusion  matrix. 

Now,  eq.  (3.15)  together  with  properties  (3.18)  may  be  modelled 

as  a  stochastic  differential  equation  if  the  sources  R  (t)  are 

a 

idealized  as  stochastic  processes,  denoted  as  R3(t),  which  have  the 
form, 


Ra(t)  =  MmnCy]  Sn(t)  (3.19) 


where 


<S  (t)  S    (t)>  =  <S  (k,t)  S  ,(k\f)> 


=  (2tt)3  6(k+k')  5(t-t»)  A  , 
  cm ' 


S  Aaa«  6(t_t,)  (3-20) 
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where  A       is  a  constant  matrix  and  the  amplitudes  M    [y]  are  fixed  by 
aa  mn 

the  fluctuation-dissipation  theorem,  eq.  (3.18).    Finally,  the 
variables  y  (t)  are  replaced  by  stochastic  variables  denoted  as 
y  (t)  =  y  (k,t)  and  the  average  over  the  distribution  function 

3.  Ot 

replaced  with  that  over  the  stochastic  sources, 

^a(t)>ss  *  Vt}>R(t)  (3'21) 
since  the  variables  y  (t)  may  be  considered  to  be  functionals  of  the 

a 

stochastic  sources. 


Long  Wavelength  Fluctuations 
If  we  restrict  our  attention  to  long  wavelength  fluctuations, 
then  the  preceding  Langevin  model  can  be  explicitly  evaluated. 
Specifically,  we  only  consider  wavevectors  such  that 


*Amfp«1.  (3.22) 


In  this  case,  the  various  functionals  appearing  in  eq.  (3.15)  may  be 

evaluated  in  terms  of  a  power  series  in  k£  .  .    If  we  truncate  the 

mfp 

expansions  at  second  order,  then  the  model  may  be  given  explicitly,5 
in  terms  of  the  usual  microscopic  conservation  laws, 


ft  9  +  1  '  1  =  0 


30 


ft 6  + 1 '  a  =  0 


It  ja  +  jk  fca6  "  °  (3'23) 


where  sa  and  tag  are  the  heat  flux  vector  and  stress  tensor, 
respectively.    These  in  turn  each  have  a  convective  part,  a 
dissipative  part  and  a  stochastic  part,  where  the  latter  represents 
the  stochastic  source, 


g  R 

s    =  eu    +  t  -u.  +  X  - —  T  +  s 
a         a       aB  6       o  3r  a 

-a 


U  -  V»B  +  P6ae  "  Vctfn  U  }  "  *o6a6  ^  +  &  (3-24) 


where  we  have  introduced  the  fluctuating  velocity,  u    =  p~1 1  the 

a  °a 

fluctuating  pressure  and  temperature,  p  and  T,  and  transport 
coefficients,  nQ»  <Q  and  X  .    These  fluctuating  thermodynamic 
functions  and  transport  coefficients  are  the  same  functions  of  y  as 
one  has  in  equilibrium;  i.e., 


p[y]  =  p  [y] 


"  PeCp.eJ  (3.25) 
where  p    is  the  equilibrium  pressure  functions  and  the  last  line 
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indicates  that  pg  depends  only  on  the  (local)  density  and  internal 
energy  density,  e,  where  e  is  given  by 


e  -  e  -  -  p     j  .  (3.26) 


The  dissipative  parts  of  the  heat  flux  vector  and  stress  tensor  are 
seen  to  be  given  by  a  fluctuating  version  of  Fourier's  law  and 
Newton's  viscosity  law,  respectively,  where 


Finally,  as  previously  mentioned,  the  stochastic  sources  are  modelled 
as  Gaussian  white  noise  with  covariances  determined  by  the 
fluctuation-dissipation  theorem,  eq.  (3.18), 


S*  =  (kRA  T2)1/2  s  (r,t) 
a         Bo  a  — 


fcae  =  (kBTTlo)1/2  W^'t}  +  (kBTKo)1  /2  (3'28) 


<s  >  =  <s.     >  =  <s„    >  =  0 
a  lag  2ag 


<3a(r1,t1)se(r2,t2)>  =  2  6a6  Sir^)  Sit^) 


<WVVS2liV(£2»t2)>  =  2  6a6  6yv  6(W  6(VV  (3'29) 


with  all  other  cumulants  equal  to  zero.    Equations  (3.23)  -  (3.29) 
completely  specify  the  fluctuating  hydrodynamics  model  which  may  be 
summarized  as  a  Langevin  equation  for  the  conserved  densities, 
y  (t),  with  a  deterministic  part  given  by  the  Navier-Stokes  equations, 
local  equilibrium  transport  coefficients  and  susceptibilities 
supplemented  by  the  stochastic  sources,  eqs.  (3.28)  and  (3.29),  which 
model  the  excluded  degrees  of  freedom.    Finally,  we  note  that  this 
model  is  given  in  the  language  of  eq.  (3.15)  in  Appendix  A  where  it 
should  be  noted  that 


—  DaJy;t]  -  o  (3.30) 

8yb 


as  is  shown  by  Zubarev  and  Morozov. 


Small  Length  Scale  Fluctuations 
In  a  dense  fluid,  the  mean  free  path  may  actually  be  smaller  than 
the  typical  atomic  length  scale,  to  be  denoted  here  as  a  (for  a  gas  of 
hard  spheres,  o  would  be  the  hard  sphere  diameter).    In  this  case,  one 
may  study  atomic  length  scale  fluctuations  with  the  present  methods. 
However,  the  evaluation  of  eq.  (3.15)  requires  care  for  while  we  still 
expand  in  (k£mfp)  we  must  keep  all  orders  in  (ka).  Quantities 
dependent  on  (ka)  carry  information  about  the  small  scale  structure  of 
the  fluid.    We  can  evaluate  (3.15)  while  retaining  the  desired 
information  regarding  the  small  scale  structure  by  taking  its  short 
time  limit  to  obtain 


33 


ft  ya(t)  -  vJyU)]  =  Ra(t).  (3.3D 


The  drift  vector  is  defined  to  be 


vaCy(t)]  =  -  (<6(y-y)  £  ya>ss/<5(y-y)>ss}|y=;(t).  (3.32) 


In  general,  this  is  a  complicated  nonlinear  functional  of  y(t). 
However,  as  is  explained  in  Chapter  5,  we  will  only  be  interested  in 
this  functional  expanded  to  second  order  in  the  deviations  from  the 
stationary  state.    Thus,  we  expand  as 


v  Cy(t)]  =  v  [y  (k)]  +  6y.(t)  •  [  -J        v  [y(t)]]    „  .  ♦  etc. 

3    °       •        b  «yb(t)    3  y0(^ 


(3.33) 


The  first  order  term  is 

sr — v  Cy(t)]  =  (<<5y<Sy>~l)ho  <5y„£y  >   .  (3.3*0 

,5y  (t)    a  ss  be      Jc    a  ss 

This  model  is  explicitly  evaluated  in  Chapter  6  and  analyzed  in 
Chapter  7.    While  it  may  seem  strange  to  study  the  short  time  limit  of 
the  generalized  Langevin  equation,  it  should  be  noted  that  similar 
short  time  equations  have  been  found  to  provide  qualitatively  accurate 
descriptions  of  dense  fluids  even  out  of  the  short  time  regime  (an 
excellent  example  is  the  Enskog  equation*^) • 


CHAPTER  4 
LINEARIZED  HYDRODYNAMIC  MODEL 

In  this  chapter,  we  linearize  the  Navier-Stokes-Langevin 
equations  described  in  Chapter  3  and  use  them  to  calculate  various 
statistical  properties  of  the  system.    We  will  do  this  in  such  a 
manner  that  the  results  apply  for  shear  rates  up  to  and  including  the 
large  shear  regime  as  previously  defined.    Of  particular  interest  will 
be  the  hydrodynamic  propagators  and  the  two  time  correlation 
functions.    We  will  discuss  the  applications  of  our  results  to  light 
scattering  experiments  and  will  make  contact  with  the  small  shear 
results  of  other  authors.    The  chapter  concludes  with  a  description  of 
various  formal  statistical  properties  of  the  system  highlighting  the 
difference  of  this  model  from  local  equilibrium. 

Linearized  Model 
The  complete  hydrodynamic  model  was  given  in  Chapter  3,  eqs. 
(3.23)  -  (3.29),  and  will  not  be  repeated  here.    The  variables  we 
choose  to  work  with  are  the  mass  density,  total  energy  density  and 
momentum.  So, 

ya  *  (p.e.j)a  (4.1) 


3^ 
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and  we  denote  the  steady  state  averages  of  these  quantities  with 
subscript  naught, 


y0a~  (p0'VVa  {H'2) 


and  the  deviations  from  the  steady  state,  the  fluctuations,  are 
denoted  by  Za, 


Za(-*t)  "  ya(-'fc)  "  y0a(-,o)  (4'3) 

where  we  also  use  internal  energy  density  and  velocity  in  place  of 
total  energy  density  and  momentum  density.    Finally,  we  scale  and 
Fourier  transform  Za(r,t),  to  obtain  the  variable  Z(k,t)  in  terms  of 
which  we  wish  to  work, 


Za(k,t)  =  /d3r  e1^  Za(r,t) 


(4.H) 


Za(k,t)  «*•  Cc^pCk.t),  c26e(k,t),  e(a)(k).u(k,t)]  (4.5) 


-2  8p0  2  -1  8pn 


where  pQ  =  PLECp0»e0]  is  the  steady  state,  macroscopic  pressure; 
h0  =  e0  +  p0  is  the  enthalPy  density  and  {e(a)  (k ) }  are  a  set  of 
pairwise  orthonormal  vectors  with  e(1)(k)  =  k/k  =  k  and  e  is  the 
internal  energy  density  and  u  =  p"1j  is  the  velocity.    The  scaling, 
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(4.6),  will  simplify  the  form  of  the  equations  below  and  we  emphasize 
that  c-j  and  C2  are  equilibrium  thermodynamic  susceptibilities. 

In  the  course  of  Fourier  transforming  and  linearizing  the 
equations,  we  find  that  we  must  transform  the  convective  terms  which 
include  the  operator  Pq1^  •  V  which  on  an  arbitrary  function  f(rO 
looks  like 

P~0\  '  Vf(r)  =  a.  .r.  JL  f(r).  (4.7) 
The  transform  of  which  is 


(4.8) 


Thus,  it  turns  out  that  the  form  of  the  linearized  equations  is 

ft  V*.t)  "  a.  .k.         Za(k,t)  +  Ha6(k;a)  Zg(k,t)  =  R^k.t)  (4.9) 

3 

where  Ra  represents  the  stochastic  sources  which,  at  linear  order,  do 
not  depend  on  the  fluctuations,  Z0.    The  details  of  the  linearization 
including  the  explicit  forms  of         and  Ra  are  given  in  Appendix  A. 
Equation  (4.9)  is  the  primary  result  of  this  section.    In  the 
following  sections,  we  will  analyze  the  hydrodynamic  (i.e.,  normal) 
modes  described  by  this  equation.    For  zero  shear,  H      is  the  usual 


matrix  of  linearized  hydrodynamics  the  eigenvalues  of  which  describe 
sound  propagation,  heat  diffusion  and  momentum  diffusion. 

An  important  feature  of  eq.  (4.9)  is  the  coupling  of  variables  of 
different  wavevectors  implicit  in  the  gradient  term.    This  unusual 
feature  has  important  consequences  on  the  evolution  of  fluctuations 
and  is  known  as  linear  mode  coupling. 

Green's  Function 
To  solve  equation  (4.9),  we  introduce  the  Green's  function 
Gag(k.,k.' ;t)  defined  to  satisfy 

ft  W^';t)  -  aijki  W~.  W^':t)  +  Haa(^a)Ga6(^,;t)  *  0 

with    G    (k,k';0)  =  (2tt)3  6(k  -  k')  5  0.  (4.10) 

In  equilibrium,  i.e.,  a  =  0,  the  solution  of  (4.10)  requires  nothing 

more  than  the  diagonal ization  of  the  hydrodynamic  matrix  H    .  Before 

ag 

solving  (4.10)  for  nonzero  shear,  we  recall  the  equilibrium 
solution. 39 

Let  the  eigenvectors  of  Hag(k_,a=0)  be  denoted  as 

l£(l)(k)  }i=1  >5  and  the  corresponding  eigenvalues  be  {A^QO}^  5  so 
that 

Ha0(k.O)  S^Oc)  =  X(i)(k)  ^i}(k).  (4.11) 
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Furthermore,  as  Hag  is  in  general  nonhermitian,  we  introduce  the 

biorthogonal  set  of  vectors,  [y,     (k ) }         ,  defined  to  satisfy 

i=i ,  b 


n(l)(k)  •  £(J)(k)  =  6    .  (4.12) 


In  terms  of  these  functions,  the  Green's  function  can  be  immediately 
determined  to  be 

5 

GaeQi.i<';t)  =    !    ^l)(k)ngi)(k')(21r)35(k-k')exp  -  [A  .  (k)  (t-tQ)] 


(4.13) 


and  t0  is  a  constant  which  can  be  fixed  by  an  additional  boundary 
condition.    The  solution  to  eq.  (4.9)  is  then  given  by 

dk' 

V*'fc)  =  J  7TT3  Ga6(-'-'  ;t)  V*'»V 


(2ir) 

dk' 
(2tt) 


dk ' 

+  It    dT  /  — T  G    (k,k' ;t-T)  R.(k',T)  (4.14) 


where  the  wavevector  integrations  extend  to  infinity. 

Because  the  hydrodynamic  matrix  Hag  is  the  result  of  a  truncation 
at  second  order  of  an  expansion  in  (k4mf  ),  the  eigenvalues  can  be 
determined  by  perturbation  theory  in  this  same  small  parameter.  In 
equilibrium,  one  obtains  two  sound  (propagating)  modes,  a  heat  mode 
and  two  degenerate  shear  modes, 


39 


*2(k) 


2 


(4.15) 


where  explicit  forms  for  the  speed  of  sound,  cQ,  the  sound  damping 

constant,  rQ,  and  the  thermal  diffusivity,  DQT,  will  be  given  later. 

For  the  present,  we  note  only  that  the  modes  all  decay  like 

exp-( constant) k  t  and  that  these  results  can  be  used  to  calculate  the 

intensity  of  scattered  light  from  an  equilibrium  fluid  (as  described 

below).    Upon  doing  so,  one  obtains  the  classic  Landau-Placzek 
44 

formula     which  is  perhaps  the  most  successful  application  of 
fluctuating  hydrodynamics.    In  fact,  comparison  of  this  result  and 
experiment  provides  one  of  the  best  known  means  of  measuring  the 
equilibrium  transport  coefficients. 

Returning  to  the  nonequilibrium  problem,  eq.  (4.10),  we  introduce 
a  shear  rate  dependent  generalization  of  eqs.  (4.11)  -  (4.12)  as 
follows: 
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'aijki  W~.  ^^a)  +  Haa(*5a)4*)(*}a)  ■  ^(A)(«Li*)ciA)(|ci«) 

n(l)(k;a)  •  i(J)(k;a)  =  6    .  (4.16) 


The  presence  of  the  gradient  operator  in  (4.16)  greatly  complicates 
the  problem.    However,  as  discussed  in  Chapter  2,  we  restrict  the  size 
of  the  shear  rate  by  requiring  that  the  rate  of  convection  be  no 
larger  than  the  rate  of  dissipation  yielding  the  condition, 


a  <  vQk2.  (4.17) 


This  choice  is  what  allowed  us  to  neglect  heating  in  these  equations, 
since  the  rate  of  heating  goes  like  a2  -  k4  (see  eq.  2.3),  and  makes 
possible  a  perturbative  solution  of  eq.  (4.16).    The  details  of  this 
calculation  are  given  in  Appendix  B.    The  eigenvalues  so  obtained  are 


A    =  -ic  k  +  -  (r  k2  +  ak  k  /k2) 
l  0        2      0  x  y 


A  =  ic  k  +  -  (r  k2  +  ak  k  /k2) 
^         U       d      0  x  y 


X3  '  D0T  k' 


2  2 
An  -  vnk    -  ak  k  /k 
4       0  x  y 


A5  -  vQk2.  (i4-18) 
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We  see  that  the  degeneracy  of  the  shear  modes  has  been  lifted  while 
the  sound  damping  constant  has  been  shifted  and  the  heat  mode  remains 
unaffected  at  this  order. 

If  we  now  write  the  Green's  function  as 


5 

Oag(k,k';t)  =    I    n^i)(k;a)5gi)(k';a)G(i)(l<,l<';t),  (4.19) 
then  we  find  that  the  function  satisfies  the  following  equation: 


Clt  '  aiiki  9k~  +  *U)  GU)(k,k';t)  =  0 

J 


Gv  '(k,k';0)  =  (2ir)J  S(k-k').  (4.20) 


This  scalar  equation  is  easily  solved  and  if  we  introduce  a  time 
dependent  wavevector  as 

k.(t)  =  k.  -  k.aj.t  =  k.   -  6.y  kxat,  (4.21) 
then  the  solution  of  (4.20)  is 


G(i)(k,k';t)  =  (2tt)3  5(k-k'(t))E(1)(k,t) 


where  E(l)(k,t)  =  exp  -  \\    dtA.(k(-T)).  (4.22) 

z0  1 

The  explicit  form  of  the  propagators  in  (4.22)  is 
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E(1)(k,t)  =  [k/k(-t)]1/2  exp{-ic0ka(t)  -  |  rQk2B(t)} 


(2)  (1  )* 

E^;(k,t)  =  EU;  (k,t) 


E(3)(k,t)  =  exp{-DQT  k28(t)} 


Ev  ;(k,t)  =  [k(-t)/k]  exp{-vQk  6(-t)} 


E(5)(k,t)  =  exp{-v0k2g(-t)}  (4.23) 


where    a(t)  =  (2ak  k)'^  {(k  (-t)k(-t)  -  k  k) 

x         y  y 


k(-t)  +  k  (-t)  sgn(ak  ) 

+  k    ln[  k  ?  /  sgn(ak  )    X  ]1 

y  x 


B(t)  =  (akxk2)"1   |k2(ky(-t)  -  ky)  +1  (k3(-t)  -  k3)} 


,2      ,  2  ,2 

k    -  kx  +  \>  (4.24) 


These  functions  behave  for  long  times  as 


o(t)  ♦  (ak  /2k)t2 
x 


8(t)  ■*  ]-  (akx/k)2t3.  (4.25) 


Because  lim  f5(t)  >  0,  the  modes  are  all  stable  as  expected  (see  Drazin 

and  Reid     and  Case     for  further  discussion  of  the  stability  of 

macroscopic  shear  flow). 

In  order  to  illustrate  the  physical  significance  of  some  of  the 

shear  rate  dependent  effects  described  by  these  equations,  consider 

the  evolution  of  a  purely  hydrodynamic  sound  excitation  along  the 

direction  of  flow.    Specifically,  we  take  Z  (k,0)  =  £      (k)  6(k-k„) 

a  —  a      —  — -0 

with  k^  =  kQx.    The  evolution  of  hydrodynamic  excitations  is  given  by 
eq.  (4.14)  when  conditionally  averaged  over  the  stochastic  force  with 
fixed  initial  condition  Z^Ck.O).    Denoting  this  average  by  a  double 
bracket  we  have  that 


dk' 

«Z  (k,t)»  =  /   -  G  „(k,k<;t)  ZQ(k',0) 

(2w)J    op  3  - 


SOS-k^t))  ^1)(k)  E(1)(k,t).  (4.26) 


Inverting  the  Fourier  transforms  yields, 

(1)  2  ?  1/U  -rnkn0(t,-a)/2 

«Za(r,t)»  =  £y;(kQ(t))(1  +  a2t2)17^  e    0  0 

x  cos{kQ[x(t)  -  cQa(t,-a)]}  (4.27) 

where  a(t,a)  and  8(t,a)  were  defined  in  (4.25)  and  x(t)  =  x  +  ayt. 
The  time  dependent  coordinate  x(t)  is  due  to  the  expected  Doppler 
shift  in  the  sound  wave  which  travels  in  a  moving  medium.  This 
kinematic  effect  may  be  suppressed  by  restricting  attention  to  the 
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y  =  0  plane  and  without  loss  of  generality,  we  may  chose  x  =  0  as 

well.    Thus,  the  remaining  shear  effects  are  due  to  gradients  of  the 

velocity  field  and  not  the  magnitude  of  the  velocity.    Figure  (4.1) 

shows  <<Za(r=0,t)>>/<<Z  (0,0)»  for  the  choice  of  parameters 
2-1 

^C0k0/r0k0^  =  10    '  whiGh  is  typical,  and  the  two  cases 
2 

(a/rQk  )  =  0  and  1.    Note  that  in  the  former  case, 

-r  k2t/2 

«Za(0,t)»/«Za(0,0)»  =  e  cos(c0k0t)  (4.28) 


which  is  the  expression  for  ordinary  sound  propagation.    Two  effects 
are  easily  seen.    First,  the  effective  sound  velocity  increases  with 
time  relative  to  that  for  zero  shear  rate  on  this  time  scale. 
Specifically, 


^.i(nM2)"!.il,[(,.A"!M]|  (..29) 
2 

where  t  =  r0kQt.  Second  is  that  the  attenuation  is  enhanced,  relative 
to  that  for  zero  shear,  as  time  increases, 

r0k0e(t'-a)      ,       1  22 

 2          =  1  +  T  a  t  .  (4.30) 

W 

These  relative  magnitudes  of  convection  and  dissipation  would  be  very 
difficult  to  reach  with  simple  fluids  but  are  possible  for  laboratory 
colloidal  suspensions. 
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Correlation  Functions 


The  correlation  functions  of  the  deviations  of  the  fluctuations 
from  the  steady  state  are  useful  for  a  number  of  reasons.    From  them 
one  may  calculate  the  dynamic  structure  factor,  which  describes  the 
intensity  of  light  scattered  from  the  system  and  transport  properties 
such  as  the  shear  rate  dependent  renormalized  viscosities  and  they 
characterize  the  deviation  of  the  system  from  local  equilibrium.  In 
this  section,  we  calculate  the  two  point  correlation  functions. 

The  two  point  correlation  functions  are  defined  to  be 


where  the  brackets  denote  the  nonequilibrium  average  which  in  the 
present  case  corresponds  to  an  average  over  the  stochastic  forces. 
Time  reversal  invariance  and  stationarity  require  that 


where  PQ  =  ±1  is  the  parity  of  Za  under  time  reversal.  Therefore, 
may  consider  t  _>  0  without  loss  of  generality. 
If  we  define 


CaB(k,fk'ft';a)  =  <Za(k,t)  Z  (k'.f )> 


(4.3D 


c„ ntk.tjk'.t' ;a)  =  C  o(k,t-t';k,0;a) 


Caf}(k.t;k',t';a)  =  P^C^Ck'  ,t;k,f  -a) 


(4.32) 


F(i)(k,t)  =  n(i)(k;a)  R  (k,t), 
—  a      -         a  — 


(4.33) 
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then,  using  eq.  (4.14),  we  find  that 


Ca6(k,t;k\0;a)  =  C    (k,k»  ;t) 


I  /"  dT/;dxM^i^l<)^J)(k')E(i)(!<,x)E(J)(J<.,T.) 
i  t  J 


l  <F(i)(k,(-T),t-T)F(j)(k'(-T'),-T')>  (4.34) 


where  we  have  taken  tQ  +  -«  in  eq.  (4.14)  and  used  the  boundary 
condition  that  Mk,-00)  is  finite.    It  is  clear  that 


<F(l)(k,t)F(J)(k',f)>  -  (2Tr)36(k+k')6(t-f)F(iJ)(k)  (4.35) 


where  the  form  of  F(lJ)(k)  is  easily  determined.  Thus, 


Ca(}(k.t;k')  =  (2lr)35(k+k'(t))  I        }  (k )      )  (-k . )  E(  1 }  (k,  t)C  0(-k';a). 

(4.36) 


The  functions  Cag(k;a)  are  related  to  the  equal  time  condition 
functions  because,  in  particular, 

caeQ<.t=0;k'  ;a)  =  (27r)35(k+k')Cae(k;a).  (4.37) 
The  explicit  form  of  these  functions  is 
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V!<;a)  =  i;  dx  I  {^i)(k)^J)(-k)E(i)(l<,x)E(J)(-k,x) 
i.j 


x  F(iJ)(k(-x))}.  (4.38) 


The  nonzero  elements  of  C  0(k;a)  are  given  in  Tables  4.1  and  4.2. 

Sore  immediate  information  may  be  gained  from  these 
expressions.    In  particular,  eq.  (4.36)  shows  that  the  two  body 
correlations  decay  in  the  same  way,  i.e.,  with  the  same  propagator 
E  *^(ktt),  as  the  macroscopic  modes.    This  represents  a  generalization 
of  Onsager's  regression  hypothesis  to  nonequilibrium  states. 
Furthermore,  eq.  (4.38)  is  a  fluctuation-dissipation  theorem  for  the 
nonequilibrium  system.    These  statements  are  more  obvious  when  cast  in 
the  form  of  differential  equations: 


Clt  "  aijki  g^WH'H'it)  +  Haa(k;a)  Coe(k,k',t)  =  0  (4.39) 
C_6aa  aijki  air  +  Haa(k-;a)]  Ca6(k-;a)  +  [V_k-;a)]  Caa(k-'a)  =  V*J 


(4.40) 


where  Rag(k_)  is  the  amplitude  of  the  fluctuating  forces, 


<R  (k,t)  R  (k»,f)>  =  6(k+k»)  6(t-f )  R  0(k).  (4.41) 

0t  p 

Equation  (4.39)  shows  more  clearly  the  linear  regression  law  while 
(4.40)  has  the  usual  form  of  a  fluctuation-dissipation  relation. 
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Table  4.1     Equal  time  correlations. 


Cpp(*;a)  =kB0  0XT[1  +  A^k.a)] 


V*:a)  =kBToVoXT[1   -n^T0  +  V*a)] 


C£e(k;a)  =  kBT0xT[x;Vv  +  [h0  "  T0  ^  +  *W±'M 


C33(k;a)  =  kBTQPQ  [1   +  YA^ (k;a)] 


C41|(k;a)  =  kBTQPQ1[1   -  A2(k;a)] 


CJ45(k;a)  =  I<bT0Pq1  A3(k;a) 


C55(k;a)  =  kBToPQ1[1  -  A4(k;a)] 
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Table  4.2    Nonequilibrium  contributions  to  equal  time  correlations, 
-ir-    kkvkv(-t)  -r0k2e(-t) 

A  (k;a)  =  aY  1 J  dt  — ^          e  ° 

,      ,         f-     kxky(_t)  -2V2S(-t) 
A2(k;a)  =  2aJ0dt    e  0 

k 

A3(k-;a)  "  a/odthrfer  F(k-'t}  + 

f-  kykv(_t)  k,  -2vnk26(-t) 

A^kja)  =  2aJ0dt  F(k,t)[  \  y(.t)  F(k,t)  ♦  jj£]e  ° 

F(k,t)  =  A(k(-t))  -  A(k) 

k  k  k 

A(k)  =   1  tan"1(-I) 

k  k  k 
x 

k(t)  -  k  -  k  •  a  t 

k    =  k    +  k 
x  z 
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Physical  Properties 

We  now  discuss  two  physical  properties  which  may  be  determined 
from  the  correlation  functions.    They  are  the  long  range  nature  of  the 
correlations  and  the  light  scattering  spectrum. 

The  density-density  equal  time  correlation  function  in  real  space 
is  easily  seen  to  have  the  form 

dk 

CnJr'a)  m  \   a  e  "  -  C  (k;a) 

=  Jl~3kBT0p^xT[6(r/£)  +  %/r  F(r/4)]  (4.42) 
where  the  isothermal  compressibility,  xT  =  P~1(|^)T,  and  where 

T  9p  T 

1  /2 

I  =  (rQ/a)        is  a  typical  length  scale  associated  with  the  shear 

rate.    The  limit  F(r/J,)=0  yields  the  usual  equilibrium  result  where 

the  6-function  is  an  idealization  of  the  atomic  scale  correlations 

which  actually  exist  in  the  fluid.    The  part  porportional  to  r/fl,  has 

been  obtained  before  via  an  expansion  to  first  order  in  the  shear 
20 

rate.       Because  of  the  form  of  F,  this  actually  requires  that 
r/l  «  1  and  yields 

F(0)  =  -a    r.r./8iTYar2  (4.43) 

*  J     ^  J 

where  Y  =  cp/cy  is  the  ratio  of  specific  heats.    This  does  not,  as  has 
20 

been  claimed,      however,  allow  us  to  conclude  that  the  correlations 
are  long-ranged  because  it  is  the  product  of  a  short  distance 
expansion.    To  determine  the  true  nature  of  the  correlations,  at  large 
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distances,  we  consider  a  somewhat  simpler  problem  than  the  large-r 
expansion  of  F(r,l).    Specifically,  we  consider 

Jdy/dz  Cpp(r)  =  r1kBT0p^xTC6(x/Jl)  +  YH(x/i)].  (4.44) 

The  part  of  the  correlation  function  due  to  the  deviation  from 
equilibrium,  H(x/4),  is  given  by 

HCx/4)  =  (2A)"1  rQ  dtf(l  ♦  tV3/2  t"1/2  (1  ♦  1  t2)"1/2 

x  exp[-(x/Jl)2(4t(1  +jt2))"1]  (4.45) 
and  for  large  (x/4)  behaves  as 

H(x/4)  *  0.828(x/A)'5/3.  (4.46) 

At  short  distances,  this  function  is  Gaussian  and  the  complete 
function  is  shown  in  Fig.  4.2.    It  is  thus  clear  that  the  correlations 
are  in  fact  long  ranged,  decaying  as  a  power  law. 

The  second  physical  property  we  wish  to  study  with  the 
correlation  functions  is  the  intensity  of  light  scattered  from  the 
fluid.    It  can  be  shown39  that  the  intensity  of  light  at  the  point  r 
with  frequency  a  due  to  the  initial  scattering  of  a  coherent  beam  of 
wavevector       and  frequency  u0  is  given  by 


bo 
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I(r,w)  -  J-  IQ  a)1,(a/2TrCQr)sin  *  S(kjQ)  (4.47) 

where  IQ  is  the  spectral  intensity  of  the  incoming  beam,  a  is  the 
polarizability  of  the  particles  composing  the  medium, 

•       s  Eq  cos<}>,  Eq  is  the  amplitude  of  the  coherent  beam, 
j£  ■  kg  ~  (uQ/CQ)r  and     =  oj  -  wQ  and  where  the  structure  factor, 
S(k_,Q),  is  defined  as 

dk  dk 

S(k,0)  -    J-dt  J   T-  e(k-k  )e(k+k?)Cn  (k1t;k„,0;a)eiQt.  (4.48) 

Thus,  we  see  that  the  density-density  correlation  function  determines 
the  intensity  of  scattered  light.    The  functions,  6(k),  are  Fourier 
transforms  of  form  factors  which  limit  the  volume  of  integration  to 
that  irradiated.    In  Appendix  C,  it  is  shown  that  S(k,Q)  can  be 
written  in  terms  of  nonequilibrium  Brillouin  and  Rayleigh  peaks: 

S(k;Q)  =  SR(k;fl)  +  SB(k;fl)  +  SB(k;-Q)  (4.49) 

where 


t&iQ)  "  2Re{[-io  +  Xl(k;a)  -  \  ^"Vft1  >n<1  \  (k;a) } 


SR(k;«)  =  2Re{[-iQ  +  A3(k;a)  -  i 


2  aijki  3kjJ  S  ci  \  cal(k;a)j. 


(4.50) 
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In  equilibrium,  the  gradients  are  not  present  and  these 
expressions  reduce  to  simple  Lorentzians.    Expansion  of  eqs.  ( 4 . 19)  to 
first  order  in  the  shear  rate  yields  previously  obtained  results; 
however,  the  general  evaluation  of  (4.19)  is  difficult.    A  quantity 
which  can  be  calculated  directly  is  the  Landau-Placzek  ratio,  R, 
defined  as 

R(k;a)  *  Jdw  SR(k ,u) /Jdw[SB(k ,M)  +  S  (k,-u)] 


(4.51) 

where  a  =  PQo^/hQc2.  In  particular,  the  integrated  Rayleigh  intensity 
is  unchanged  from  equilibrium, 

Jdw  SR(k,u>;a)  =  ( Y  -  1)/Y  (4.52) 

where  Y  =  ratio  of  specific  heats.    An  explicit  form  of  (4.51), 
obtained  using  the  expressions  for  the  equal  time  correlation 
functions  given  in  Appendix  C,  is 


R(k,a)  =  (Y  -  1 )/[1  +  YA  (k;a)] 


(4.53) 
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where  A    =  [C    (k;a)  -  C    (k;0)]/C    (k,0).    In  Fig.  4.3,  we  plot 
i         pp  -  pp  —  pp  - 

AR  =  [fHk;a)  -jUk  ;0)  ]  /R(k  ;0 )  as  a  function  of  9  for 

2 

k  -  kQ(cosex  +  siney)  and    a/rQk0  =  1.    The  small  shear  rate  expansion 

of  these  results  is  in  agreement  with  earlier  work  by  Kirkpatrick  et 
20 

al.       The  result,  eq.  (4.53),  is  unexpected  in  that  it  states  that 
the  integrated  intensity  of  the  Rayleigh  peak  is  unaffected  by  the 
shear. 

Probability  Densities 

In  this  section,  we  recast  our  previous  results  in  the  language 
of  the  probability  distributions  of  the  fluctuations.    In  particular, 
the  probability  distributions  are  more  useful  than  the  description  in 
terms  of  fluctuations  for  the  calculation  of  higher  order  correlation 
functions  or  for  averages  of  more  general  functions. 

The  probability  and  joint  probability  density  are  defined  by 

P(Z,t)  =  <6(Z(t)  -  Z)> 


P(Z,t;Z\t')  =  <6(Z(t)  -  Z)  6(Z'(f)  -  Z')>  (4.54) 

where  the  averages  are  taken  over  the  stochastic  forces,  {Ra},  and 
over  a  specified  ensemble  of  initial  values,  {z(0)}.    The  6-functions 
are  short  for 


<5(Z(t)  -  Z)  =    n    6(Z  (k,t)  -  Z  (k)). 

.        a  —  a  — 

a,k 


(4.55) 


We  begin  the  evaluation  of  eq.  (4.51)  by  calculating  a  related 
quantity— the  conditional  probability  density  defined  by 


W(Z,t;Z\t')  =  P(Z,t;Z\f)/P(Z',f).  (4.56) 

Because  the  system  is  being  modelled  by  a  Markovian  process,  the 
conditional  distribution  is  independent  of  the  initial  values, 
lza(0)}         because  of  the  neglect  of  heating,  the  system  is 
stationary  so  W(Z,t;Z\t')  =  W(Z,t-t'  ;Z » ,0) .    So,  without  loss  of 
generality,  we  consider 

W(Z,t';Z',0)  =  <6(Z(t|Z»)  -  Z>R  (4.57) 

where,  as  indicated,  the  average  is  only  over  the  stochastic  forces 
and    Z(0|Z')  =  Z'.    From  eq.  (4.14), 

za(t|Z')  =  G(t)Z«  +  /J  d-r  G(t-x)  R(T)  (4.58) 

where  we  have  suppressed  internal  summations.    An  integral 
representation  of  the  6-function  in  (4.57)  and  use  of  (4.58)  yields 

W(Z,t;Z',0)  =  JdA  eU,CZ-G(t)Z,]<exp[-iA.  /Jdx  G(t-t)  R(t)> 

=  JdA  expjiACZ  -  G(t)Z']  -  A2M} 
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=  ldet[TTM(t)]}"1/2exp{-  |[Z-G(t)Z']  •  M~1(t)  •  [Z-G(t)Z-]} 

(4.59) 

where  the  second  line  follows  from  the  first  upon  use  of  the  Gaussian 
statistics  of  the  source  term.    The  matrix  M(t)  is  given  by 

M(t)  =  /Jdi/Jdr'  G(t-t)  •  <R(x)R(x')>  •  GT(t-T').  (4.60) 

Equations  (4.39)  and  (4.40)  (the  regression  law  and  fluctuation- 
dissipation  theorem)  allow  us  to  re-express  (4.60)  in  terms  of  the 

correlation  function  matrix,  C(t)  =  C  D(k,t;k' ,0) , 

otp  —  — 


M(t)  =  C(0)  -  C(t)   •  C_1(0)   •  CT(t).  (4.61) 


We  can  now  immediately  write  down  the  probability  and  joint 
probability  densities, 

P(Z,t)  =  JdZ«  W(Z,t;Z',0)  P(Z\0) 

P(Z,t'Z',t')  =  W(Z,t;Z',f)  P(Z'.t'),  (4.62) 

where  P(Z',0)  is  the  initial  probability  density.  From  (4.61),  we  can 
determine  the  stationary  state, 
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P    (Z)  =  lim  P(Z,t) 
ss 


{det[TrC(0)]}"1/2  exp[-  |  Z  •  c"1(0)  •  Z]  (4.63) 


which  follows  from  the  limits, 


lim  C(t)  =  0  =>  lim  M(t)  =  C(0)  (4.61) 
t+°°  t+o> 

and,  lim  G(t)  =  0.    Thus,  as  expected,  because  the  driving  process  was 
t+<*> 

Gaussian  and  the  equations  linear,  the  distribution  of  fluctuations  is 
also  Gaussian.  Furthermore,  eq.  (4.63)  shows  that  the  same  asymptotic 
stationary  state  is  reached  regardless  of  the  initial  condition. 

All  of  these  results  could  have  been  obtained  if  we  had  worked 
with  the  Fokker-Planck  equation  rather  than  the  fluctuation 
equation.    For  completeness,  we  write  down  the  Fokker-Planck  equation 
for  P(Z,t)  which  is  easily  obtained.    For  example,  by  differentiating 
(4.62)  with  respect  to  time, 

It  P(Z'U  =  "I    STlky  CJ0ct(Z't)  +  J1a(Z't)]  (M-65) 

OL  f  K  0t 

where  the  drift  term  is 


Joo(z-t)  =  Caijki  atT  -  Lae(!l^)ZB(k)]P(Z,t) 


(4.66) 


and  the  diffusion  term  is 


J1a(Z't}  "5  P(Z't)-  (4'67 

6 

The  matrix  R    (k),  the  diffusion  matrix,  is  just  the  correlation 

up 

matrix  of  the  stochastic  forces.    Equations  (4.65)  -  (4.67)  are  the 
usual  results  for  the  relationship  between  a  linear  Langevin  equation 
and  the  corresponding  Fokker-Planck  equation. 


CHAPTER  5 

NONLINEAR  FLUCTUATING  HYDRODYNAMICS 

In  this  chapter,  we  will  consider  the  consequences  of  the 
nonlinear  terms  occurring  in  the  general  fluctuating  hydrodynamics 
model.    We  begin  by  considering  a  naive  form  of  renormalization  which 
illustrates  in  an  intuitive  way  the  concept  of  renormalization.  Using 
this  method,  we  will  derive  the  lowest  order  renormalization  of  the 
shear  viscosity  explicitly  deriving  the  nonanalytic  dependence  on 
shear  rate.    In  the  following  sections,  we  will  present  a  more  formal 
development  of  renormalization  which  will  result  in  a  self-consistent 
formulation  of  the  nonlinear  problem.    The  chapter  concludes  with  the 
connection  between  the  naive  and  formal  renormalization. 

Naive  Renormalization 
Consider  a  simple,  one  dimensional,  nonlinear,  stochastic 
equation  of  the  form 

It  X  +  L0X  +  V3  =  f(t)  (5.1) 
with  a  Gaussian,  white  stochastic  source  f(t)  having  correlations, 

<f(t1)f(t2)>  =  f06(trt2)  (5.2) 
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and  all  higher  order  cumulants  equal  to  zero.    Furthermore,  let  the 
average  of  x(t),  <x(t)>,  be  denoted  by  x0(t).    Then,  the  equation  of 
motion  for  xQ(t)  is  found  by  averaging  eq.  (5.1), 

ft  X0(t)  +  L0X0(t)  +  N0<x3(t)>  =  0.  (5.3) 


If  we  write  x(t)  =  x  (t)  +  6x(t),  then  (5.3)  becomes 


ft  X0(t)  +  CL0  +  3N0<(6x(t))2>]xQ(t)  +  NQx3(t)  =  -  NQ<(6x(t))3>. 

(5.4) 


We  see  that  the  equation  for  xQ(t)  has  the  same  form  as  that  for  the 
fluctuating  variable,  x(t),  except  that  the  "bare"  coefficient  LQ  has 
been  renormalized  by  the  nonlinearity  to  become  the  "dressed"  or 
renormalized  coefficient  L, 


L  =  LQ  +  3NQ<(6x(t))2>.  (5.5) 


If  we  subtract  eq.  (5.4)  from  (5.1)  and  solve  the  resulting  equation 
for  6x  neglecting  the  nonlinear  term,  and  use  this  to  compute  L,  we 
find  that 


L  =  L0  +  3N0f0/2L0  +  o(N0}  (5.6) 


where  we  have  imposed  stationarity ,  |^  <x(t)x(t)>  =  0,  to  fix  the 
equal  time  correlation  functions. 
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The  previous  calculation  can  be  reformulated  in  terms  of 
correlation  functions.    The  equation  of  motion  for  the  two-point 
function,  <6x(t)6x(0)>  =  c(t),  is  found  from  (5.1)  to  be 

|^  C(t)  +  LQC(t)  +  N0<(6x(t))36x(0)>  =  0  (5.7) 

where  we  have  used  the  condition  <f(t)<5x(0)>  =  0.    To  lowest,  NQ  =  0, 
order,  we  find  that  the  Laplace  transform  of  C(t),    C(z),  is 


C(z)  =  C(0)/(-z  +  Lq)  +  o(NQ)  (5.8) 


which  has  a  pole  at  the  bare  transport  coefficient  LQ.    If  we 
calculate  the  correlation  function  in  (5.7)  to  lowest  order,  we  find 
that 

<(6x(t))35x(0)>  -  J^dx^/0.^  e'Vt-V^-^-Vt-VV* 


<f(T1  )f  (T2)f(T    )f(Tj|)>    +   0(N  ) 


3  ft  j     fO  j        ^V^V  'L0(t-T2)  2 

3  )-J^).J^2  e  e  2  fQ 2  +  o(NQ) 


1       1       "L0t  2 
32L^2L^e        f0  +0(V 

3f 

2L^  c(°)(t)  +  o(V"  (5.9) 
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So, 

N  f 

|^  C(t)  +  (LQ  +  \  -—■)  C(t)  =  o(NQ2)  (5.10) 

and 

C(z)  =  C(0)/(-z  +  LQ  +  |  NQf0/L0)  +  o(NQ2)  (5.11) 

which  has  a  pole  at  the  "renormalized"  coefficient  defined  in  (5.5). 
This  coincidence  is  easy  to  understand.    Suppose  that  at  time  t-0,  the 
distribution  function,  p(x),  is  perturbed  from  the  steady  state 

function,  p  ,  in  the  following  way 
s 

p(t-0)  =  pg  e'h6x.  (5.12) 
The  auxiliary  field,  h,  might  be  chosen  such  that 

<x(0)>h  *  <x(0)>ss  =  xQ,  i.e.,  to  yield  some  particular  initial  value 

for  the  macroscopic  variable.    For  small  deviations  from  the  steady 

state  (i.e.,  h  small  so  <x(0)>    -  <x(0)>      is  small),  we  find  that 

n  s  s 

<x(t)>  -  xQ  =  C(t)h  +  o(h2) 

=  (C(t)/C(0))[<x(0)>h  -  xQ]  +  o(h2).  (5.13) 


,  the  correlation  function  acts  as  the  Green's  function  for  small 
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deviations  of  the  macroscopic  variables.    They  must,  therefore,  share 
the  same  linear  relaxation  times.    A  general  proof  of  this  property 
for  Markovian  processes  has  been  given  by  Dufty  et  al. 

In  the  same  way,  we  can  derive  an  expression  for  the  renormalized 
shear  viscosity.    If  we  take  the  average  of  eqs.  (2.45),  we  require 
that  we  obtain  the  macroscopic  hydrodynamic  equations.    For  shear 
flow,  the  only  nontrivial  equation  is  the  heat  equation, 


where  nR  is  the  renormalized  viscosity.  Comparing  (5.14)  with  the 
average  of  eqs.  (2.45)  yields 


where  we  have  used  temporal  and  spatial  translational  invariance  to 
set  overall  derivatives  of  the  correlation  functions  to  zero.  The 
third  line  of  (5.15)  follows  directly  from  the  momentum  equation. 

Using  our  results  from  the  previous  chapter,  eq.  (5.15)  can  be 
written,  to  lowest  order  in  the  nonlinearity ,  in  terms  of  the  equal 


(5.14) 


-2  3 
nR  =  ~a      <tij(ll»t)  7^-  u.(r,t)> 

J 

=  a'2  <Ui(r,t)  ^-  t.j(r,t)> 
J 

=  a'2  <u.(r,t)  ~-  Cp(r,t)  u  (r,t)  u.(r,t)]> 
j  J 

-2  3 

=  a      <[3r~  ui(H't)]  P(£.t)  u^r.t)  u  (r,t)> 


(5.15) 
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time  correlation  functions, 


nR  =  a  1pQ<ux(r,t)  u  (r,t)> 


=  -(2ir)      Jd  kCkxkyc33(k;a)  +  k^c^k  ;a)  -  k^Ckia)].  (5.16) 


This  equation  for  the  lowest  order  renormalization  of  the  viscosity  is 
an  example  of  mode  coupling.    That  is,  the  original  nonlinearities  in 
the  equations  of  motion  have  coupled  two  different  hydrodynamic  modes, 
analogous  to  phonons  in  solid  state  physics,  and  the  correction  to  the 
viscosity  is  due  to  the  interaction  of  these  modes,  analogous  to  the 
phonon-phonon  scattering  contribution  to  the  coefficient  of  thermal 
diffusion  in  solids.    The  coupling  of  the  modes  in  (5.16)  is  buried  in 
the  equal  time  correlation  functions.    It  can  be  seen  by  referring  to 
eq.  ( -4 . 38) .    Contributions  also  exist  due  to  three  mode  coupling,  four 
mode  coupling,  etc. 

Upon  evaluation,  we  find  that 


1  /2 

\  =  nQ  +  r^a  (5.17) 


with 


0  v  0 

Thus,  it  is  seen  that  a  nonanalytic  dependence  of  the  viscosity  on 
shear  rate  arises  as  a  result  of  mode  coupling.    This  is  in  contrast 
to  the  analytic  dependence  one  would  obtain  from  the  standard  Chapman- 
Enskog  solution  of  the  Boltzmann  equation  (or  the  Enskog  equation;  see 
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Chapter  6).    This  result  agrees  with  that  of  Ernst  et  al.16  and  with 
that  of  Yamada  and  Kawasaki16  (although  their  numerical  evaluations 
are  incorrect).    Thus,  this  result  indicates  the  equivalence,  with  the 
present  assumptions,  of  this  model  and  the  others  appearing  in  the 
literature. 

Formal  Renormalization 
In  this  section,  we  wish  to  develop  a  more  general  method  of 
renormalization  than  the  naive  scheme,  allowing  the  calculation  of  all 
renormalized  properties.    To  this  end,  we  begin  by  formulating  the 
nonlinear  model  in  a  more  general  notation, 


It  ym  +  vm[y]  +  DmnCy]FnCy]  =  MmnCy]5n  (5'19) 


where  vm[z]  represents  the  nonlinear  Euler  terms  in  eq.  (3.23)  (i.e., 
the  terms  V  •  j  =  V  ■  (pu)  in  (3.23)  and  ev    and  v  j    in  (3.210).  The 

Ot  Qt  3 

matrix  D^Cy]  is  called  the  diffusion  matrix  and  F  [y]  are  the 
thermodynamic  forces  which  are  the  variables  thermodynamically 
conjugate  to  the  ym's.    The  product  DmnFn  represents  the  dissipative 
terms  in  (3.24).    The  reason  for  returning  this  form  of  the  equations 
is  that  in  their  derivation  of  fluctuating  hydrodynamics,  Zubarev  and 
Morozov    show  that  thermodynamic  consistency  imposes  the  fluctuation- 
dissipation  theorem, 


(5.20) 
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"thermodynamic  consistency"  meaning  that  the  equilibrium  distribution 
function  is  a  solution  to  (5.17). 

Equation  (5.17)  is  highly  nonlinear  making  a  general  analysis 
difficult.    However,  in  equilibrium  it  has  been  found  that  truncating 
eq.  (5.17)  at  second  order  in  the  nonl inearities  provides  a  qualita- 
tively correct  description  of  such  nonlinear  phenomena  as  long  time 
tails  (i.e.,  the  algebraic  decay  of  equilibrium  correlation  func- 
tions).   In  fact,  in  equilibrium,  the  dominant  mode  coupling  effects 
arise  from  the  Euler  terms  in  (5.19).    We  shall,  therefore,  consider 
only  these  nonlinearities  in  the  rest  of  the  discussion.    It  can  be 
shown  that,  to  order  (ak),  this  is  equivalent  to  taking  both  Dmn  and 
Mmn  t0  lowest  order  (i.e.,  Djy]  -  Dmn[yQ])  thus  preserving  the 
fluctuation-dissipation  theorem,  eq.  (5.20),  in  this  truncated  model 
of  eq.  (5.19).    Furthermore,  Fn[y]  is  expanded  to  first  order  in  the 
fluctuations  in  order  to  preserve  the  usual  linearized  transport  laws. 

With  these  approximations,  and  expanding  the  variables  ym  about 
their  steady  state  averages  yQm,  eq.  (5.19)  takes  the  form, 


3t  za(k,t)  ♦  Lab(k,t)  zb(k,t)  +  eVab0<k.k1.k2)(zb(k1ft)  zc(k2,t) 


-  <zb(k1,t)zc(k2,t)>) 


-  R  (k,t) 


(5.21) 


where  L 


ab 


+  Hah^'a)  is  the  linear  operator  studied  in 
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the  previous  chapter  and       ~  y&  -         (details  are  given  in  Appendix 
B).    The  vertices,  which  arise  solely  from  the  Euler  terms,  have  the 
form 

=  *abo(*»v  5(*-w  (5-22) 

where  the  explicit  form  of  A^k.gJ  is  given  in  Appendix  B.  The 
constant  e  is  introduced  as  a  notational  convenience  in  the 
development  of  a  perturbative  expansion  in  the  nonlinearity . 

Equation  (5.21)  can  be  rewritten  using  the  Green's  function 
introduced  in  Chapter  4  as 

«a(k,t)  -  Z^W)  ♦  «jSdT{0Mf<k.Vt-T)  Vatbc(qrq2,q3) 

X  (zb(VT)  Zc(%'T)  "  <2b(a2'0)zc(%'0:)>)J  (5.23) 
where  a  sum  over  repeated  wavevectors  is  implicit  and  where 

Za°)(*;t)  "  Gaa'(^;t)  za'(4>°>  +  /Jdx  O^.C^t-t)^, (£.,). 

(5.24) 

Using  (5.23),  a  perturbative  solution  to  eq.  (5.21)  is  easily 
constructed. 

From  eq.  (5.21),  we  can  immediately  construct  the  equations  of 
motion  of  the  two-point  correlation  functions, 
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C  .  (k,k';t)  =  <z  (k,t)z(k'0)>, 
ao   a  —       b  — 


ft  Cab(fc,K'.t)  +  Laa,(k;a)  Ca,b(k,k';t)  +  M^Oc.k' ;t)  =  0  (5.25) 


where  we  have  defined 

Mab(k,k';t)  =  V^Ck.q^q^  <zc(qi,t)zd(q2,t)zb(k',0)>  (5.26) 

and  have  used  the  fact  that  <R  (k,t)  z  fk',t')>  =  0  for  t*<t 

a  —        b  — 

(causality—the  source  at  time  t  cannot  affect  a  fluctuation  occurring 

at  an  earlier  time).    The  process  of  renormalization  occurs  when  we 

define  a  self-energy  (memory  function)  T  w(k,k':t)  as 

ab  

Ma>b(k,k';t)  =  /JdT  Iaa,(k,ki;t~r)  C^U^  ,k'  ;x) .  (5.27) 

The  significance  of  this  is  that  eq.  (5.24)  becomes  a  linear  equation 
for  the  correlation  functions  and  effective  viscosities,  etc.  as  would 
be  measured  in  the  laboratory  may  be  identified. 

Laplace  transforming  (5.24)  in  conjunction  with  (5.27),  one 
obtains 

Z  5ab(^';w)  +  Laa'(^;a)  £a'b(^'  ;w)  +  <a'(^rw)  -H' 

"  cab(H.k';t=0)  (5.28) 
and  the  Fourier  transformed  self-energy  may  be  computed  from  (5.27), 
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Eab(-,k-1 JM)  =  \c(^'^;w)  C~  ob^S.k'iw).  (5.29) 

As  discussed  in  the  previous  chapter,  C    (k,k»;t)  =  C    (k ,  t  )<5(k-k '  ( t) ) 

ab  —  ab  —        —  — 

and  it  is  easily  verified  that,  order  by  order  in  perturbation  theory, 
Mab(-'-' ;t)  S  M    (k,t)6(k-k'(t))  which  requires  that 
lab^.k'.t)  s  Iab(k,t)6(k-k'(t)).    Thus,  eq.  (5.28)  becomes 


ft  Cab(^'t}  +  D^>Cab(k,t)  +  Haa,(k;a)Ca,b(k,t)  +  eM^Oc.t)  =  0 

(5.30) 


where  we  have  introduced  D  =  -a  ^  -|-  .  The  equation  for  the  self- 
energy  becomes 


Mab(^t}  =  /odT  IM.^t-T)CaIb(k(T-t)fT)) 


which  shows  the  self-energy  to  be  an  operator  in  k-space.  Laplace 
transforming  (5.31 )  yields  an  expression  for  the  self-energy, 

J0dt  e       ^aa-(ii;t)e       -  Mab(*;w)5ba'(*;w)- 


3  ^aa'(-'w)  (5-32) 
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and  the  Laplace  transform  of  (5.30)  is 

[wD]Cab(k,w)  +  Haaf(k,a)Ca,b(k.w)  +  e/0"dtIaa,(k;(O)e-(w+D)(T)5a,b(k,w) 

=  Cab(k,t=0).  (5.33) 

Equations  (5.32)  and  (5.33)  are  the  primary  results  of  this 
section.    The  linear  mode  coupling  which  was  the  unique  feature  of  the 
linear  problem  is  seen  to  play  an  important  role  here  as  well.  The 
operator  nature  of  the  self-energy  is  a  result  of  the  linear  mode 
coupling  and  expresses  the  fact  that  new  gradient  couplings,  as  well 
as  the  expected  new  hydrodynamic  couplings,  arise  as  a  result  of  the 
process  of  renormalization.    This  complication  is,  however,  removed  in 
the  context  of  perturbation  theory  in  e  because  we  have  an  expression, 
to  lowest  order  in  e,  for  the  operator  (w+D)  acting  on  the  correlation 
functions  (c.f.  eq.  (5.33)).    It  is  seen  that  the  self-energy  function 
shifts  the  poles  of  the  correlation  function  renormalizing  the 
transport  coefficients  as  in  eq.  (5.11). 

Application  of  Formal  Renormalization 
We  now  wish  to  apply  the  renormalization  formalism  to  the 
calculation  of  the  renormalized  shear  viscosity  and,  in  so  doing,  to 
make  a  connection  with  the  naive  renormalization  method.    As  a 
simplification,  we  will  restrict  our  attention  to  the  contribution 
arising  from  the  shear  modes.     In  this  case,  the  irreversible  part  of 
the  stress  tensor,  t    ,  can,  in  general,  be  written  as 
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"  vijlm  9^  Jm  '  (5'3J4) 


where  v  m  is  the  viscosity  tensor  which,  in  the  unrenormalized 
theory,  has  the  Navier-Stoke's  form, 


Vijlm  "  ^6ii6jm  +  5im6jl  -  §  +  (5.35) 


where  v  and  k  are  the  usual  kinematic  and  bulk  viscosities, 
respectively.    The  macroscopic  equation  for  the  energy  density  is 


ft  e  +  3^  C(e  +  p)Po1j.  +  t.1p*1j1  ♦  s.]  =  0  (5.36) 


For  macroscopic  shear  flow,  using  eq.  (5.34),  the  heating  equation 
becomes 

It  6  +  Vijlm  aijalm  =  °-  (5.37) 

Thus,  the  renormalized  viscosity  calculated  earlier  from  the  heating 
equation  corresponds  to  the  yxxy  element  of  the  general  viscosity 
tensor. 

We  now  calculate  the  self-energy  to  second  order  in  the 
nonlinearity  (i.e.,  to  second  order  in  e).    Using  eqs.  (5.23)  and 
(5.26),  we  find  that 
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Mab(k,k';t)  =  eVacd(k,q1,q2)<z^0)(qrt)z^0)(q2>t)zb(k'>0)> 


+  2e2/otdT{Vacd(^a1.a2)Gce(a1,a3;t-T)Vefg(q3,q4,%) 


x  <40)(^'T)zg0)(%'T)zd0)(%'T)zb(^''o)>J +  °(e3) 


(5.38) 


where  the  factor  of  two  arises  from  the  symmetry  of  the  vertices. 
Using  eq.  (5.24)  and  causality,  this  reduces  to 

Mab(k,k';t)  -  «VMd(M1.fle)0oe,(fl1.fl3lt)Oddi(flB.a4«t) 

x  <zc,(q3>0)zd,(qljf0)zb(k',0)> 


Gff.(£l4.aj;;T)Ggg,(q5,q^;T)Gdd,(q2,q2;t) 

x  <zf(£r0)zg(35,0)zd(22,0)zb(k',0)>}  +  o(e3). 

(5.39) 

In  equilibrium,  the  equal  time  fluctuations  are  Gaussian  distributed 

making  the  o(e)  term  in  eq.  (5.39)  zero  and  allowing  easy  evaluation 
2 

of  the  0(e  )  term.    In  the  steady-state,  the  equal  time  fluctuations 
are  not  necessarily  Gaussian  distributed.    Stationarity ,  in  the  steady 
state,  allows  us  to  write 
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za(k;0)   =  (k;0)   ♦  «/!.dt{0M((k,ai»-T)Vafbo(ai.fl2.a3) 

x  (zb(q2;x)zc(q3;T)  -  <zfe (£g ; Oz^ ; T ) » } 

(5. MO) 

where 

za0)(k;0)  =  J^dx  Gaa,(k,a;-T)Ra(ai-T)  (5.41) 

and 

za0)(k,-t)  =  J~*dx  Gaa,(k,q;-t-T)Ra(qrx)  (5.12) 

where  in  (5.42),  t>0.    Thus,  the  equal  time  fluctations  may  be 
determined  perturbatively  and  are  Gaussian  distributed  to  zerotn  order 
in  e. 

Equations  (5.39)  -  (5.42)  provide  us  with  enough  information  to 
calculate  the  self-energy  function  from  eq.  (5.28).    The  result  is 
that 

£aa.Q<.i<";w)  -  e/0-dt  e-Wt{Vacd(k,2l,a2)Gcc,(ai,a3;t)Gdd,(a2,ai|;t) 

*  <2c.(%.0)zd,(24,0)zb(k',0)>}c^a;"1(l<',k«;z) 
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P  r  oo  —  wfc  r 

+  Wodte  {Vacd(^ra2)Vega,(q3,q1(,!<") 

x  Gde(a2'a3;t)Gcd'(^'^;t)}<Zd,(q^0)z  (q^,0)> 


+  o(e3).  (5.H3) 


To  recover  our  previously  calculated  lowest  order  mode  coupling 

contribution  to  the  viscosity,  we  first  replace  the  equal-time 

correlations  in  (5. 43)  with  their  local  equilibrium  values.    The  term 

proportional  to  e  in  (5.43)  then  vanishes  and  (5.43)  becomes 

k  T 

Eaat(k,k«;w)  --f^Jo"  e"WtlAacd(l<,q)Aeyat(k(-t)-q(-t),-q(-t)) 

-D(k)t 

x  ^e^'S't^Qd.^.'tJe  }6(k+k"). 

(5.44) 

The  momentum-momentum  self-energy  upon  evaluation  is  found  to 
have  the  form 


y-*z)  -  kikj  w-*z)  (5-_*« 


which  identifies  it  as  the  renormalization  of  the  viscosity  tensor. 
If  we  write  v.  =  v<°>m  ♦  Sv'..^,  then  the  long  wavelength,  long 
time  hydrodynamic  viscosity  appearing  in  eq.  (5.37)  is 
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«vww  =  lim  I   (k,w)  (5.16) 

k+0 
w+0 


yxxy     k+Q  "yxpxpy  - 


from  which  the  "naive"  result  is  recovered. 

The  replacement  of  the  static  correlation  functions  by  their 
local  equilibrium  values  requires  some  comment.    For  a  stationary 
process,  the  correlation  of  the  stochastic  sources  is  related  through 
fluctuation-dissipation  theorems  to  the  static  correlations.    The  mode 
coupling  contributions  to  the  static  correlations  may  thus  be 
interpreted  as  a  renormalization  of  the  amplitudes  of  the  stochastic 
sources.    By  neglecting  the  equal  time  three-point  function  in  (5.42), 
we  were  left  with  the  propagator  renormalization. 

There  is  another  way  to  understand  the  replacement  of  the  static 
correlations  in  the  present  calculation  by  their  local  equilibrium 
values  in  order  to  find  agreement  with  the  "naive"  calculation.  The 
only  information  about  the  "statics"  that  was  used  in  the  "naive" 
calculation  was  the  stochastic  source  correlations.    But,  as  explained 
in  Chapter  3,  these  are  approximated,  in  the  unrenormalized  theory,  by 
their  local  equilibrium  values.    It  is,  therefore,  reasonable  that  the 
same  information  must  be  used  in  both  calculations  to  obtain 
agreement.    We  speculate  that  if  the  first  order  renormalized  forces 
were  used  in  the  naive  calculation,  one  would  obtain  the  same  result 
as  is  obtained  here  without  the  local  equilibrium  replacement. 


CHAPTER  6 
SHORT  TIME  MODEL 


In  this  chapter,  we  discuss  the  explicit  form  of  the  short  time 
model  of  hydrodynamic  fluctuations  introduced  in  Chapter  3.    We  recall 
that  this  model  is  intended  to  describe  hydrodynamic  fluctuations  on 
length  scales  as  small  as  atomic  lengths.    We  begin  by  discussing 
general  properties  of  the  distribution  functions.    This  is  followed  by 
a  derivation  of  explicit  expressions  for  the  necessary  matrix  elements 
in  terms  of  the  reduced  distribution  functions.    The  final  two 
sections  describe  the  approximations  introduced  to  evaluate  these 
matrix  elements  and  the  explicit  form  of  the  resulting  equations 
governing  the  time  evolution  of  hydrodynamic  fluctuations. 

General  Properties  of  the  Reduced  Distribution  Functions 
In  this  section,  we  introduce  two  properties  of  the  distribution 
functions  which  will  be  of  use  in  determining  the  general  nature  of 
the  expressions,  eqs.  (3. 31 ) -(3. 34) ,  to  be  evaluated  below  and  which 
will  aid  in  their  evaluation.    These  properties  are  the  stationarity , 
in  time,  and  the  modified  translational  invariance  of  the  distribution 
functions. 

As  discussed  in  Chapter  2,  an  external  nonconservati ve  force  is 
applied  in  computer  simulations  of  sheared  fluids  to  counterbalance 
viscous  heating  of  the  system.    The  significance  of  this  is  that  it 
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yields  a  steady  state.    If  we  denote  by  p    (r)  the  distribution 

ss 

function  (density  matrix)  of  the  steady  state  system,  depending  on  all 
of  the  phase  space  variables  r  =  lSLi»£ilJ.i«  tnen  we  can  write  the 
steady  state  Liouville  equation  as 


0  "  ft  pss(r)  +  £"-P(r) 


=  «.p(r)  (6.2) 

where  £_  is  the  Liouville  operator.    The  explicit  form  of  £  is 
generally  given  by 


-      n=1m£n      9an      2^n,         9qn  3£n  +  £ext  (6'3) 

where  V(|qJ1-qJl,|)  iS  the  pair  potential.    For  hard  spheres,  the 
potential  is  singular  rendering  eq.  (6.3)  ill  defined.    As  discussed 
in  Chapter  3,  however,  eq.  (6.3)  can  be  replaced  by  a  well  defined 
expression  which  produces  a  statistically  equivalent  dynamical  system 

-         N  1 

£~  =    h  £n  *  3q~  +  2    ^    ^_(n,n')  +  T_(n',n)]  ♦  f  (6.4) 
n=l  -^n        n*n*  ext 

where  we  have  included  a  term  representing  the  externally  applied 
thermostat  alluded  to  above.    The  explicit  form  of  the  T  operator  and 
its  various  adjoints  and  time-reversed  forms  are  given  in  the 
literature  and  will  not  be  repeated  here.    The  point  we  wish  to  make 
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is  that  (6.2)  allows  the  £+  operator  to  be  moved  in  averages,  i.e., 
for  any  two  phase  functions,  A(r)  and  B(r),  it  can  be  shown  that 

<A(r)(£   B(r))>     =  <(£  a(D)  B(r)>   .  (6.5) 

ss  -  ss 

The  subscript  on  £_  refers  to  the  fact  that  T+  goes  to  its  time 
reversed  form  T_. 

Another  important  property  of  the  distribution  functions  is  their 
modified  translational  invariance.    Clearly,  because  the  local  flow 
velocity  varies  spatially,  the  distribution  function  cannot  be 
translationally  invariant.    That  is,  if 

-  Ei^i  "  R  (6.6) 
for  an  arbitrary  displacement  R  ,  then 


However,  if  we  make  the  transformation 

*  * 
Ei-S-i  -  S.L  ~  maijRj.  lt  -  R  ,  (6.8) 

then  we  expect  the  distribution  function  to  be  invariant  since  the 
macroscopic  state  is  characterized  by  the  local  velocity  and  the  local 
velocity  gradient  and,  while  (6.6)  does  not  preserve  these  properties, 
(6.8)  does.    Mathematically,  it  is  easy  to  verify  that  (6.8)  is 


83 


consistent  with  the  Lees-Edwards  boundary  conditions  while  (6.6)  is 
not.    We,  therefore,  look  for  solutions  of  Liouville's  equation  which 
are  functions  of  the  variables 


a1<3  ~=  £i  -  aj  (6.9) 


which  are  invariant  under  (6.8), 


£1'%      *     £i*'a*j-  (6.10) 


Henceforth,  we  will  write  p    (r)  a  p    (fD'.  a  }) 

ss  Kss  l^i'  %  j  J 


Explicit  Form  of  Matrix  Elements 
The  short  time  model  is  given  in  Chapter  3  as 

atV^t)  +  La6(k,k';a)6£6(k',t)  =0  (6.11) 

where 

V^';a)  ■  <(£+ya(k)>*yY(!<")>33  (g~\6,k«,k'  (6-12) 


and 
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gYR(k",k')  =  <5yv(k")6yQ(k')>    .  (6.13) 

IP  T  p  —  SS 


The  variables,  5y  (k ) ,  are  again  taken  to  be  the  fluctuations  of  the 
conserved  densities  around  the  steady  state: 


fiy,  (k)  =  P(k)  =  —  I  me1-'^ 
/V  n 

6y2(k)  =  6e'(k)  -  —  I  -  |  k  T  )e^^n 

2  /V  n  2    B  o 

«y2+6(!<)  =  6j:(k)  =  4  I  (p:)r  cs.14) 

/V  n     n  3 


which  are  summarized  as 


6ya(k)  =  4  I  J  (j.)  e^-Sn,  (6.15) 
✓V  n 


where  V  is  the  volume  of  the  system. 

We  now  define  the  reduced  distribution  functions, 


nnf(n)(X1 Xn}  "  (MhnTT  H+1         dXn  pss(r,)  (6'l6) 


where       «♦  (an»£n)  and  we  note  that  the  Jacobian  of  the  change  of 
variables  Xr     x;  is  unity.    Using  (6.16)  and  the  independence 
of  pgs  on  the  center  of  mass,  we  find  that 


BoB<!i»!S,>  "  (2lr)3V_16(k+k')  gae(k),  (6.17) 


where 

gag(k)  -  n  /dp'  f(1)(£.)  ;a(£.)  y£.) 

The  explicit  form  of  £  y  (k)  is  easily  found  to  be 

+  a  — 

£+C  (k)  -  i  4  I  H  •  ^  i  (R')e^^n 
a  /v  n         m    a  n 

+  ?4    I    [T  (n,n*)  +  T  (n',n)]  Y  (p')e^- 

"  ki3ij  3T  Ja(^  +  <extV^ 
3 

From  this,  it  follows  that 

<(£+y  (k))6y  (k')>      =  (2TT)3«5(k+k')M  Q(k) 
a  p         ss  —  —     ag  — 

"  kiaij  9kT 

+  ^xtV^V^^ss' 

where 
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Ma6(k)  =  lnk./dE}^io(E[)iB(Ef)f(1>(E{) 


+  \  n2/d£'d^d£12(CT+(1,2)+T+(2J)])ia(E|)yE.)f(2)(E.>£.;ai2) 

♦  °2  /dEidE2dai2(CT+(1'2)+T+(2'1)])^(Ei)V£2)el~  1112 
x  f(2)(£1',E';a12) 


+  §7  JdEj  d^d^dq,  2dq1  3  ( CT+  ( 1 ,  2)  +T+  ( 2 , 1 )  ] )  Ya (£• ) ) y 6      )  e1"  £l  3 

(3) 

x  f    (a,  .£2 '£3 ;a12-a13)«  (6.21) 


Approximations 

The  evaluation  of  eqs.  (6.21)  and  (6.18)  require  knowledge  of  the 
one-  and  two-body  reduced  distribution  functions  for  the  steady 
state.    In  general,  these  functions  cannot  be  determined  explicitly 
and,  so,  we  introduce  two  important  approximations. 

The  first  of  these  is  the  neglect  of  momentum  correlations,  which 
is  to  say  we  assume  the  two-body  distribution  to  have  the  form 


f(2)(X',X')  "=  f(2)(£l',E2;£l2) 


=  ^"(P^f^Cp^gCajg).  (6.22) 
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where  f  1^(j>' )  is  the  one-body  reduced  distribution  function,  which 
cannot  depend  on       due  to  the  modified  translational  invariance,  and 
the  function  g(q.12)  is  the  steady  state  analogue  of  the  pair 
correlation  function  (i.e.,         2^dSl 2  is  the  Probability  of  finding  a 
particle,  2,  in  the  volume  dq^  2  if  there  is  a  particle,  1,  at  position 
£1).    In  equilibrium,  this  approximate  form  becomes  exact.  However, 
away  from  equilibrium,  we  expect  that  particle  velocities  are  indeed 
correlated  (i.e.,  (6.22)  implies  that  <£{£^>sa  =  0  but  in  general  this 
need  not  be  the  case).    This  approximation,  which  is  analogous  to 
Boltzmann's  Stosszahlansatz, ^  is  standard  in  the  development  of 
kinetic  equations.    In  particular,  it  is  used  to  derive  the  Enskog 
equation  which  is  known  to  give  a  qualitatively  accurate  description 
of  fluids  near  equilibrium.1*3    With  this  approximation,  the  three-body 
term  in  (6.21)  only  contributes  for  6=1.    in  Appendix  D,  it  is  shown 
that  this  term  can  be  written,  using  an  exact  identity  following  from 
stationarity,  in  terms  of  two-body  functions. 

Having  adopted  the  form  (6.22),  we  must  still  determine  the  one- 
body  distribution  and  the  pair  correlation  function.    The  second 
approximation  we  use  is  to  replace  the  steady  state  pair  correlation 
function,  g(q12),  by  the  local  equilibrium  pair  correlation  function, 
8LE(|l12|)  =  8le.(q12^      The  pair  correlation  function  is  related  to 
the  density-density  equal  time  correlation  function  and  knowledge  of 
one  of  these  allows  determination  of  the  other.    This  approximation, 
which,  like  the  previous  one,  is  used  in  the  derivation  of  various 
kinetic  equations,  may  be  thought  of  as  a  first  guess  to  be 
"renormalized"  using  the  density  autocorrelation  function  calculated 
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from  the  resulting  equations.    This  approximation  then  becomes  the 
first  iteration  of  a  self-consistent  calculation.    In  fact,  these  two 
approximations  require  that  f(1)(p_.f)  satisfy  an  Enskog  equation  which 
is  the  usual  dense  fluid  analogue  of  the  Boltzmann  equation. 

With  these  approximations,  our  model  is  completely  defined.  The 
one-body  distribution  can  be  obtained  from  the  first  equation  of  the 
BBGKY  hierarchy  (i.e.,  by  integrating  eq.  (6.2)  over  coordinates 
X£...X^  and  using  eqs.  (6.4)  and  (6.22)): 

^4"f(1)(1,)  =  /dV-^2)f(1)^)f(1)(2)gLE(ai2) 

+  sources.  (6.23) 

The  sources  in  eq.  (6.23)  arise  from  the  nonconservati ve  forces  in  the 
Liouville  equation  (eq.  6.4).    However,  it  can  be  shown  that  they 
contribute  to  second  order  in  the  shear  rate  and  so  are  neglected  in 
the  present  calculation  which  is  only  carried  out  to  first  order  in 
the  shear.    While  the  previous  approximations  served  to  construct  a 
model,  we  now  make  an  approximation  with  respect  to  solving  that 
model;  we  will  only  work  to  first  order  in  the  shear  rate.    Thus  we 
expand  f ' 1  \ 

f  "teP  =  fLE(£i)(1   +  ■•♦!<£{)  +  -..)  (6.24) 

where  fLE(£»)  is  the  local  equilibrium  distribution  function  for 
uniform  shear  flow, 
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-3  /2  2 

fLE(Rj)  ■  n(2miTkBTo)  J      exp  -  (pj c/2mkBTQ) .  (6.25) 

With  these  approximations,  one  is  able  to  determine  ^(jaj)  in  eq. 
(6.25).    As  this  result  exists  in  the  literature,  3  we  only  quote  it 
here: 

a/1V  =  -^jVoVPliPij  "  I  Pl,2<Sij) 
nk  =  x  C1  +  it  TO*x]nB 

%  =  ~T  (m/ire0)1  72  (6-26) 

1 6o 

and  x  -  SLE(a).    The  constant,  rig,  is  the  kinetic  part  of  the 
viscosity  as  calculated  from  the  Boltzmann  equation.    It  should  be 
noted  that  we  have  only  sought  the  steady-state  solution  of  eq. 
(6.22). 

With  these  approximations,  the  matrix  g    (k)  may  be  determined  to 

otp  — 

be 

g1 1 (k,a)  =  nS(k) 
g22(k,a)  =  |  n(kBTQ)2 

gct8(^a)  =  pOkBT06a6  "  mVai  i+aii }        a=i+2'  8=j+2=3,4,5 


ga3  =° 


all  others 


(6.27) 
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while  the  matrix  M  0(k,a)  can  be  written  as 
dp  — 


ik-£». 

M      =  -n<          y  (p')  y  (p')> 


+  XoB  "  Sag  (6.28) 


where  <...>E  denotes  an  average  over  fE(£*)  truncated  to  first  order 
in  a.    The  collisional  contribution,  Iag,  is 


Ia6(k,a)  =  -n2X/d£l'd£'d£12  ya(£.)  T_(12)f(1)(£pf(l)(£j.) 

x  ^B(£')  +  Y6(£')e       12  -  531e  ~H2]  (6.29) 


and  the  source  term,  S    ,  is 

a8 


Sa6  =  ^extV^V"^  ~  V'^ss^as'  (6.30) 
It  is  shown  in  Appendix  D  that         -  a2,  allowing  us  to  neglect  it 
here. 


Hydrodynamic  Equations 
Once  the  various  integrals  have  been  performed,  we  find  the 
following  form  for  the  short  time  hydrodynamic  equations,  where  the 
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stars  indicate  dimensionless  quantities  with  lengths  scales  to  a  and 
times  to  tg, 


D  p    -  ixk«j_    =  0 


*     *  #        2        *    ~  "  * 

D.6e    +  [DT(x)x    +  a  c,k  k  ]6e 
t  T  1  x  y 


-  iCxh(x)ka  -  a*(c2  *  %)(ky6ax  ♦  kx6ay)  -  a%3kxkyka]6j^ 


o(a2) 


— * 

D*6j*  -  i[xP  (x)k    -  a*  ^st-t  (k  6      +  k  6  )]6p* 
c    01  Pa  oyS(x)      y  ax       x  ay  F 


-* 

iCxP  (x)k    -  a*(c,  +  X)(k  6      +  k  6  )]6e* 

e       a  1         tt  y  ax       x  ay 


+  |  [v*(x)x26a6  +  (v'*(x)-v*(x))x2kake]6j* 


+  a  [0,(6    6„    +5    6„  )  +  c  k  k  <5 

6    ax  By       ay  flx'       7  x  y  afl 


+  cftk  (k  {.     +  k  5„  )  +  n  (k  6      +  k  6  )k 
8  a    y  Bx       x  By         9    y  ax       xay  B 

+  C10kxkykak6]6JB  +  a\x6By6JB  =  ^  <6-3D 


where  x  =  ka  and, 
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yx)  - 1  d  -  j0(x))/x2 


v'*(x)  =  I  (1  -  j.(x)  +  21  (x))/x2 


v*(x)  =  I  (1  -  J0(X)  -  j2(x))/x2 


h  (x)  =   r  (1  +  3yj1(x)/x) 

36y^  1 


P*(x)  =  -2-^  S"\x) 
P  36y^ 


Pe(x)  =  |  (1  +  3yj1(x)/x)  (6.32) 


are  wave vector  dependent  generalizations  of  the  thermal  conductivity, 
viscosity,  shear  viscosity,  enthalpy,  density  derivative  of  pressure 
and  energy  derivative  of  pressure.    The  constants  y  =  -  Tma3x 

and  v    =  2Hy  ^1  +  5         v    is  a  diroensionless  form  of  the  kinetic 
contribution  to  the  viscosity,  The  functions  jn(x)  are  nth-order 

spherical  Bessel  functions  with  j^x)  =  (sin  x)/x.  The 


functions  c.  =  c.(x;a)  are  given  in  Table  6.1.  Finally, 

*  3     .*         2  *       a         *  s 

P    =  Pa  ,  j    =  jtE/a,  etc.  while  Dfc  -        -  a^k.         la  the 

convective  derivative  in  terms  of  dimensionless  timeV  =  t/t  and 

*  E 
a    =  at£. 

At  this  point,  we  make  a  few  comments  about  the  general  nature  of 
these  equations.    The  first  of  eqs.  (6.31)  is  the  familiar  continuity 
equation  describing  the  conservation  of  mass.     The  second  two 
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equations  describe,  respectively,  the  conservation  of  energy  and 

momentum.    However,  due  to  the  fact  that  the  effects  of  atomic 

structure  are  included,  they  are  more  complicated  than  at  small 

wavevector.    These  complications  disappear  if  we  take  x  <<  1  and 

expand  to  second  order  in  x  in  which  case  we  get  a  short  time  version 

of  the  linearized  Navier-Stoke's  equations.    We  note  that  in  h,  P  p 

'    e'  p 

and  c2  we  have  picked  up  second  order  contributions  to  the  pressure  in 
the  long  wavelength  (small  x)  limit.    To  illustrate  the  meaning  of  the 
new  couplings,  we  have,  in  Table  6.2,  given  the  expansions  of  the 
coefficients  Cj  to  second  order  in  x.    It  is  seen  that  all  of  the 
collisional  contributions  to  the  nonequilibrium  couplings  are  of  third 
order  or  higher  in  the  gradients  except  for  c2  (recall  that  these 
coefficients  are  all  multiplied  by  the  shear  rate,  a,  which  is  a  first 
order  gradient).    Thus,  they  do  not  appear  in  the  ordinary  Navier- 
Stokes  equations  but  rather  in  the  higher  order,  so  called  Burnett, 
equations.    They  are  always  present  in  the  analysis  to  be  presented 
below,  however,  as  we  keep  all  orders  in  ka. 

Also,  these  equations  contain  wavevector  dependent  "transport 
coefficients,"  eq.  (6.32),  which  are  the  high  frequency  limit  of  the 
general  wavevector  and  frequency  dependent  transport  coefficients, 


* 

DT(x)  =  lim  D  (x,w),  (6.33) 


etc.  Wavevector  dependent  susceptibilities  such  as  the  pressure 
derivatives  also  appear.    An  important  feature  of  the  transport 


Table  6.1    Coefficients  appearing  in  hydrodynamic  equations. 


c1  - 

* 

C2  " 

2           1  -« 

J3) 

* 

C3  = 

-2C1  ♦  |  v*]J3 

* 

12    2  -* 

57  y  ci  ♦  4v  ](Jl 

+  J3) 

* 

C5  = 

12    2  -* 

-~  y  ci  ♦  nv  ]j3 

* 

c6  " 

2             2  -* 

-  f  yd  +  -  v  ](1 

2  J2; 

X 

°7  = 

-6y[1  +  |  v  ]  j  /x 

* 

c8  " 

* 

°7 

* 

°9  ■ 

* 

C7 

C10  " 

-6y[1  +  |  v*][j2 

7  •  i 
-x  J3] 

Table  6.2    Small  wavevector  limits  of  the  hydrodynamic  coefficients. 


w  v  20  -*  2 


3tt 


*  2  r        1  -* 

°2  *  "  T5  [1   +  2  v  ]x 


* 

+  0 


*         4      2r  -*  2 

175^  y  C1  +  4v  ]x 


cc  ■»  0 


*  1  2  -*  2 

°6  *  ~  35  yCl  +  7  v  ]x 


*  2  2  — *  2 

C7  *  "  35  y[1  +  7  v  ]x 


#  * 

c8  *  c7 


C9  *  c7 


c10  *  0 
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coefficients  is  that  they  decay  faster  than  x~2  for  large  x  which 
contributes  to  the  mode  softening  of  generalized  hydrodynamics  which 
was  alluded  to  in  the  introduction  and  which  shall  be  explicitly 
presented  in  the  next  chapter. 


CHAPTER  7 
ANALYSIS  OF  THE  SHORT  TIME  MODEL 

In  the  previous  chapter,  we  obtained  the  equations  of  a  short 
time  hydrodynamic  model  valid  at  large  wavevector.    In  this  chapter, 
we  study  the  hydrodynamic  modes  of  this  model  and  compare  them  to  the 
more  familiar  long  wavelength  hydrodynamic  modes.    The  primary 
characteristic  of  equilibrium  generalized  hydrodynamics  is  a  softening 
of  the  modes  at  large  wavevector  and  we  shall  find  that  this  is 
further  enhanced  at  large  shear  rates  to  the  point  of  creating  an 
instability.    We  shall  also  find  the  appearance  of  new  propagating 
modes  at  large  shear  rate. 

Equilibrium  Generalized  Hydrodynamics 
Before  studying  the  effects  of  shear  on  the  small  wavelength 
hydrodynamics,  we  review  equilibrium  generalized  hydrodynamics  and 
compare  it  to  ordinary  hydrodynamics.    If  the  shear  rate  is  set  equal 
to  zero,  the  hydrodynamic  matrix  calculated  in  the  previous  chapter  is 
found  to  be  a  function  of  ko,  where  k  is  the  magnitude  of  the 
wavevector  and  o  the  hard  sphere  diameter,  with  no  dependence  on  the 
direction  of  k.    it  is  also  easy  to  see  that  the  shear  mode,  momentum 
fluctuations  in  a  direction  perpendicular  to  the  wavevector,  decouple 
and  is  doubly  degenerate  as  in  the  small-k  limit  (c.f.,  eq.  (4.15)). 
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However,  while  the  small  k  modes  grow  as  k  ,  the  generalized  modes 
soften  at  large  k  and,  in  fact,  the  generalized  heat  mode  is  two 
orders  of  magnitude  smaller  at  ko  -  7  than  an  extrapolation  of  its 
small  k  quadratic  behavior  would  predict.    For  this  reason,  the 
relaxation  time  of  hydrodynamic  fluctuations  on  these  length  scales  is 
much  longer  than  one  would  expect  and  this  explains  how  these 
microscopic  modes  can  affect  macroscopic  properties.    Figure  7.1  is  a 
plot  of  the  real  parts  of  the  equilibrium  sound,  heat  and  shear  modes 
respectively  for  no3  =  0.88.    Figure  7.2  shows  the  equilibrium  sound 
velocity  as  a  function  of  (ka).    These  results  agree  well  with  the 
more  detailed  calculation  of  deSchepper  and  Cohen  and  support  the 
extension  of  the  present  model  to  nonequilibrium  situations. 

The  most  striking  feature  of  the  generalized  modes  aside  from 
their  softening  is  that  the  sound  mode  becomes  purely  dissipative  for 
5.75  <  ko  <  6.5.    In  other  words,  for  this  range  of  wavevectors,  sound 
does  not  propagate.    Finally,  we  should  note  that  these  modes  were 
calculated  using  the  Carnahan-Starling  expression  for  x  =  g(a)  and  the 
Percus-Yevick  expression  for  the  static  structure  factor,  S(k) ,  with 
an  approximate  Verlet-Weiss  correction.116 


Nonequilibrium  Generalized  Hydrodynamics 
If  we  now  vary  the  shear  rate,  we  find  several  striking  changes 
in  the  hydrodynamic  modes.    Before  discussing  these,  however,  we  must 
note  that  the  results  to  be  presented  here  were  obtained  for  the  case 
k  =  ky.    This  choice  of  wavevector  eliminates  the  gradient  term  in  eq. 
(6.20)  and  reduces  the  problem  of  determining  the  modes  to  a  simple 
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eigenvalue  problem.    The  case  of  wavevectors  not  along  the  y-axis  can 
be  studied  using  the  methods  presented  in  Chapter  4.    However,  the 
hydrodynamics  at  large  wavevector  is  complicated  even  for  the  choice 
of  wavevector  used  here  and  is  representative  of  the  effects  to  be 
found  for  other  wavevectors  so  that  the  simplifying  choice  seems 
best.    What  is  more,  the  results  of  Chapter  4  indicate  that  the 
slowest  decaying  modes  will  be  those  with  k  along  the  y-axis  and  so 
any  instability  in  the  system  should  occur  for  a  lower  shear  rate  than 
for  any  other  direction  of  wavevector.    Thus,  we  will  henceforth  take 
k  =  ky  unless  otherwise  stated. 

Figures  7.3  -  7.8  show  the  real  parts  of  the  sound,  heat  and 
shear  modes  for  a*  =  atE  =  0.1  through  0.64  for  n*  =  na3  =  0.88.  A 
peculiar  feature  of  these  modes  is  that  above  a*  =  0.05,  there  are  two 
sets  of  propagating,  or  sound,  modes  for  certain  ranges  of  ko.  One 
set  evolves  from  the  equilibrium  sound  modes  and  is  a  shear  rate 
dependent  version  of  these  modes.    The  other  pair  of  propagating  modes 
is  the  result  of  a  shear  induced  mixing  of  the  heat  mode  and  one  of 
the  shear  modes.    As  the  figures  show,  the  propagating  modes  appear 
and  disappear  and  mix  with  one  another  as  ka  is  varied. 

For  the  density  no3  -  0.88,  Fig.  7.8  shows  that  an  instability 
exists  for  ka  *  6.0  and  a*  =  0.64.    m  Fig.  7.9,  we  have  plotted  the 
critical  shear  rate  (i.e.,  the  minimum  shear  rate  for  which  'an  insta- 
bility exists)  as  a  function  of  density.    We  believe  this  instability 
to  be  associated  with  the  transition  to  the  so-called  string  phase 
observed  in  computer  simulations.    Also  shown  in  Fig.  7.9  are  the  two 
values  of  the  critical  shear  rate  as  reported  by  Erpenbeck.12 
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Fig.  7.9    Critica^  shear  rate,  ac*,  as  a  function  of  reduced  density, 
n*  =  no  .    Also  shown  are  Erpenbeck's  values  for  the 
transition  (o). 
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Numerically,  Erpenbeck's  values  are  about  ac*  =  0.4  and  0.8  for  n*  - 
0.704  and  0.88,  respectively.    (There  is  some  uncertainty  in  these 
numbers  as  the  transition  is  not  sharp.)    Our  values  for  these 
densities  are  ac*  =  0.63  and  0.85,  respectively.    Furthermore,  S. 
Hess,  in  private  conversation,  has  reported  that  the  transition  to  the 
string  phase  is  preceded  by  a  transition  from  the  amorphous  fluid 
phase  to  a  layered  phase  in  which  the  hard  spheres  are  partitioned 
into  planes  lying  parallel  to  the  x-z  plane.    This  is  consistent  with 
the  critical  wavevector  pointing  in  the  y-direction  since  this 
indicates  homogeneity  in  the  x-z  plane  with  variation  in  the  y- 
direction.    Finally,  the  eigenvector  associated  with  the  unstable  mode 
is,  in  our  calculation,  almost  a  pure  density  mode  while  the 
partitioning  into  planes  could  be  viewed  as  a  phase  transition  in 
which  the  macroscopic  density  goes  from  being  a  constant,  p(r)  -  p  , 
to  being  constant  in  the  x-z  plane  and  varying  periodically  in  the  y 
direction,  p(r)  -  p^  sin    iry/d  (where  d  is  the  interplane  spacing). 

The  instability  we  observe  is  attributable  primarily  to  a  single 
element  in  the  hydrodynamic  matrix:    the  coupling  of  momentum 
fluctuations  in  the  x-direction  into  the  equation  for  momentum 
fluctuations  in  the  y-direction  (i.e.,  a  coupling  of  transverse  to 
longitudinal  momentum  fluctuations  since  k  =  ky).    This  is  easily 
ascertained  in  that  when  this  coupling  is  set  to  zero,  by  hand,  the 
minimum  in  the  equilibrium  heat  mode  remains  unchanged  as  the  shear 
rate  is  increased.    From  eq.  (6.31 ),  the  coupling  in  question  is  seen 
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to  be  a*(c6*+c9*)  which,  from  Table  6.2,  is  found  to  be  proportional 
to  a*(ka)    in  the  long  wavelength  limit.    It  thus  is  of  third  order  in 
the  gradients  in  this  limit  (n.b.,  a  -  |VuJ)  and  does  not  occur  in 
the  Navier-Stoke's  equations  but,  rather,  in  their  generalization  to 
include  third  order  gradients  which  are  known  as  the  Burnett 
equations.    This  explains  why  the  instability  is  not  seen  at  small 
wavevectors  (c.f.,  Chapter  H  in  which  the  long  wavelength  equations 
are  shown  to  be  stable). 

There  is,  however,  a  second  instability  at  small  (although  not 
asymptotically  small)  wavevectors  as  seen  in  Figs.  7.3  -  7.9.  We 
believe  that  this  second  instability  is  not  seen  in  the  simulations 
due  to  the  fact  that  our  calculations  were  performed  for  an  infinite 
system  while  the  simulations  are  necessarily  done  for  finite 
systems.    Applicability  of  our  results  requires  that  Ak/k  <<  1 
where  Ak  -  2n/L,  where  L  is  the  length  of  the  box  (for  a  finite 
system,  the  wavevector  must  be  integral  multiples  of  the  "unit 
wavevector"  Ak).    For  the  simulations,  L3  -  N/n  where  N  is  the  number 
of  particles  and  n  the  density.    Thus,  applicability  of  our  results 
requires  that  (ka)»0.25  for  n*=0.88  and  N=4000  (the  smallest  Ak 
occurring  in  Erpenbeck's  simulations).    For  wavevectors  on  the  order 
of  Ak,  we  expect  that  the  boundary  conditions,  which  we  have  not 
explicitly  included  in  our  calculation,  become  important.    In  fact,  it 
is  easy  to  understand  physically  both  the  nature  of  the  instability 
and  the  reason  it  is  not  seen.    It  is  found  that  this  instability 
arises  from  the  inertial  (k-independent)  momentum-momentum  coupling  in 
eq.  (6.31).    However,  one  can  show  that  a  zero-wavevector  instability 
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arises  from  this  coupling  which  is  associated  with  center  of  mass 
motion.    This  zero-wavevector  instability  is  unphysical,  being 
prohibited  by  boundary  conditions,  and  we  conjecture  the  same  is  true 
of  the  small  wavevector  instability. 

A  similar  calculation  has  been  performed  by  Kirkpatrick^7  and 
Nieuwoldt  (hereafter,  KN) .    Two  primary  differences  exist  between  our 
calculation  and  that  of  KN.    The  first  is  that  KN  neglect  temperature 
fluctuations.    We  find  in  our  model  that  the  neglect  of  temperature 
fluctuations  drastically  alters  the  hydrodynamic  modes  and,  in 
particular,  gives  rise  to  an  unstable  mode  at  small  shear  rate  and 
both  large  and  small  wavevector.    Thus,  this  approximation  is 
unjustifiable  in  our  model. 

The  second  difference  between  our  calculation  and  that  of  KN  is 
that  KN  begin  with  the  Enskog  equation  rather  than  the  Liouville 
equation.    They  therefore  approximate  the  nonequilibrium  pair 
correlation  function  with  its  equilibrium  form  from  the  beginning. 
This  means  they  have  not  built  in  information  about  the  stationary 
state  and  as  a  consequence  cannot  use  stationarity  to  eliminate  the 
three-point  function  as  we  do.    They  are  thus  forced  to  use  an 
approximation  to  the  three-point  function  (specifically,  they  use  the 
Kirkwood  superposition  approximation)  which  we  do  not  need. 

The  result  of  these  approximations  is  that  KN  sees  an  instability 
but  at  second  order  in  the  shear  rate.    There  is  an  inconsistency, 
however,  in  that  KN  do  not  evaluate  the  effects  of  the  thermostat 
which  appears  explicitly,  at  second  order  in  the  shear  rate,  in  the 
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equations.  Thus,  the  present  calculation  uses  fewer  assumptions  than 
that  of  KN. 

Finally,  it  would  be  very  useful  to  obtain  some  idea  of  the 
completeness  of  the  description  of  fluctuations  presented  here.  The 
simplest  check  would  be  to  study  how  well  the  hydrodynamic  modes 
described  here  conform  to  those  generated  in  the  simulations.  For 
instance,  qualitatively,  one  might  check  whether  there  are  indeed  two 
sets  of  propagating  modes  at  intermediate  shear  rates. 


CHAPTER  8 
CONCLUSIONS 


The  object  of  this  thesis  has  been  a  description  of  fluctuations 
and  transport  in  a  sheared  fluid.    As  was  previously  discussed,  this 
system  is  particularly  amenable  to  theoretical  studies,  both 
analytical  and  simulational ,  but  is  difficult  to  study 
experimentally.    Experimental  studies  of  more  complicated  systems  are 
possible  and  many  of  the  phenomena  seen  in  the  simulations  of  sheared 
simple  fluids  are  also  seen  in  laboratory  experiments  on  these  more 
complex  systems.    This  agreement  lends  credence  to  the  simulations 
which  we  thus  use  in  place  of  experiment  as  a  source  of  empirical 
data. 

We  began  our  theoretical  study  by  restricting  our  attention  to 
long  wavelength  fluctuations.    These  were  described  using  fluctuating 
hydrodynamics  generalized  to  the  nonequil ibrium  state.    The  linear 
problem  was  solved  under  less  restrictive  conditions  than  had 
previously  been  used  which  allowed  us  to  include  states  far  from 
equilibrium  in  our  analysis.    Using  our  knowledge  of  the  hydrodynamic 
modes,  we  investigated  various  statistical  and  physical  properties  of 
the  fluid  including  the  behavior  of  the  time  dependent  correlation 
functions,  the  static  correlation  functions,  the  light  scattering 
spectrum,  and  the  distribution  function.    The  most  important 
consequences  of  the  shear  were  due  to  the  appearance  of  a  gradient  in 
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wavevector  space  in  the  linear  equations  (the  linear  mode  coupling 
term).    This  term  brings  about  a  much  faster  decay  of  fluctuations  and 
correlation  functions  than  in  equilibrium  (an  e       decay  as  opposed  to 
e  t).    We  also  find  long  ranged  correlations  in  the  fluid. 
Specifically,  the  equal  time  one  dimensional  correlation  of  densities 
decays  as  x"5/3  indicating  that  the  radial  distribution  function 
decays  as  r"11/3  for  iarge  r. 

We  then  preceded  to  analyze  the  effect  of  the  nonlinearities. 
These  manifest  themselves  through  a  renormalization  of  the 
microscopic,  bare,  transport  coefficients  and  susceptibilities.  The 
renormalized  shear  viscosity  was  calculated  using  a  "naive"  method  and 
found  to  be  in  agreement  with  previous  calculations.    This  important 
concurrence  indicates  the  equivalence  of  the  presence  model  with  those 
studied  previously.    A  formal  discussion  of  renormalization  in  this 
model  was  also  presented.    The  formal  scheme  allows  us  to  calculate 
renormalized  properties  other  than  the  viscosity  as  well  as  providing 
a  systematic  means  of  studying  higher  order  contributions. 

In  the  final  two  sections,  we  studied  the  derivation  and  analysis 
of  a  model  of  small  wavelength  phenomena.    It  was  found  that  the 
hydrodynamic  modes  were  drastically  modified  by  the  shear  and,  in 
particular,  that  a  large  shear  instability  was  found.    The  instability 
was  proposed  as  the  mechanism  behind  the  transition  to  the  string 
phase  observed  in  simulations. 

There  is  still  much  work  to  be  done  on  this  problem.  In  the 
fluctuating  hydrodynamics  problem,  the  nonlinear  problem  is  still 
relatively  unexplored.    At  a  practical  level,  the  complexity  of  the 
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correlation  functions  as  determined  from  the  linearized  model  makes 
the  calculation  of  the  renormalizations,  even  to  lowest  order,  a 
difficult  numerical  problem.    It  would  be  of  particular  interest  to 
calculate  the  renormalized  shear  viscosity  from  the  formal  theory 
without  using  local  equilibrium  correlation  functions  (c.f.,  the  end 
of  Chapter  5).    More  generally,  the  shear  induces  non-Gaussian 
statistics  for  the  static  fluctuations  making  this  problem 
qualitatively  different  from  most  of  those  previously  studied.  The 
consequences  of  this  difference  remain  unexplored. 

In  the  short  time  model,  several  questions  remain.    One  of  the 
most  important  concerns  the  shear  rate  expansion.    It  is  found  that 
the  shear  rate  is  scaled  to  two  different  times:    the  Enskog  time,  tE, 
and  the  Boltzmann  time,  tQ  =  xt£,  where  for  a  dense  fluid  x  -  5.  The 
highest  shear  rate  considered  is  atg  -  0.6  (near  the  instability)  or 
atE  -  0.1.    Thus,  the  expansion  in  atE  is  quite  reasonable  while  that 
in  atB  is  questionable.    However,  the  atB  expansion  arises  solely  from 
the  expansion  of  the  local  equilibrium  distribution  function.    If  this 
expansion  were  not  made,  the  elements  of  the  hydrodynamic  matrix  would 
have  to  be  calculated  numerically.    Such  a  calculation  would  provide  a 
more  quantitative  description  of  the  instability. 

It  would  also  be  useful  to  understand  the  importance  of  momentum 
correlations.    When  we  assumed  a  product  form  for  the  distribution 
function  (c.f.  eq.  6.22),  we  implicitly  neglected  momentum 
correlations.    In  equilibrium,  this  leads  to  a  50$  error  in  the 
viscosity  as  calculated  from  the  Enskog  equation.    We  expect  that  an 
inclusion  of  these  effects  would  produce  quantitative  changes  in  our 
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results.    However,  because  the  instability  is  due  to  the  appearance  of 
a  new  coupling,  not  present  at  long  wavelengths,  one  might  expect  the 
qualitative  result  of  the  existence  of  an  instability  to  be  robust 
while  the  quantitative  results  as  to  the  critical  shear  rate  must  be 
expected  to  be  sensitive  to  changes  in  the  model. 

Finally,  the  methods  used  to  study  the  long  wavelength  model 
could  also  be  applied  to  the  short  time  model.    One  might  then 
calculate  correlation  functions,  renormalized  transport  coefficients, 
etc.    These  calculations  can  be  expected  to  produce  results  which  may 
be  compared  to  data  compiled  from  simulations. 


APPENDIX  A 
LINEARIZATION  OF  THE  LANGEVIN  EQUATION 

The  fluctuating  hydrodynamics  model  was  detailed  in  eqs.  (3.23)  - 
(3.29).  We  begin  the  linearization  of  this  model  by  casting  it  in  the 
form  given  in  eq.  (5.12), 

Tt\  +  +  Dmn[y]  Fn[^]  =  MmnCy]  h'  U'1) 

This  is  accomplished  by  identifying, 

V  [y]  =  V  j 
p  a  a 

Ve[y]  =  Va(eja/p) 

Vy]  "  Vp\&  +  WP)'  (A. 2) 

FpCy]  =  -  (uCy]  -  \  j2/p)/T[y] 

Fe[y]  =  1/TCy] 

F^Cy]  -  -  ja/(PT[y])  (A.3) 
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and  all  elements  of  Dmn  zero  except 


Dee  "  VM{T[y]ACy]«ai|  +  Wylj^/pjv 


D.      -  D  , 
J  e  ej 
v  v 


In  these  equations,  T[y]  =  Tg[y]  and  y[y]  =  MgCy]  where  TgCy]  is  the 
equilibrium  temperature  functional,  etc.    The  y-dependent  transport 
coefficients  are  similarly  defined  in  terms  of  their  equilibrium 
analogs.    Finally,  in  (A.4)  we  have  abbreviated 

naBuv  =  nAaBlJV  +  <  «a6  «yv  (A.5) 
With  Aa6yv  33  in  eq*  (3-^).    The  matrix  Mmn  is  zero  except 

Me,1a6  "  '  VJ«  T^U2  VP 
Me,2a6  =  "  V«W  T^U2  VP 
Me,3a  =  " 
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and 


51ctB  =  S1a8 


C2aB  =  S2aB 


C0     ■  S  .  (a  7) 

3a       a  v.*.  f ; 


From  eq.  (3.29),  we  find  A      =0  except 

mn 


1aB,1uv  aBuv 


A2a8,2yv  =  Vuv 


A3a,3B  =  V  <A-8) 

With  these  definitions,  it  is  easy  to  verify  that  the  fluctuation- 
dissipation  relation,  eq.  (5.13),  is  verified. 

Although  this  model  could  be  studied  as  it  stands,  we  wish  to 
consider  a  restricted  model  in  which  only  quadratic  nonlinear! ties 
arising  from  the  convective  terms,  i.e.,  those  in  (A. 2),  are  kept. 
Thus,  the  restricted  model  is  obtained  by  making  the  substitutions; 


1 


D  [y]  ♦  D  [y_] 
mn  J         mn  J0 


6F 

FmCy]  +  F[y  ]  +  —21  |  6y 
m  m    0  jy  'y 


Mmn[y]  *  MmnCyO]  <A-9> 

where  Sy  =  y  -  yQ  and  the  convective  terms  to  second  order  (except  the 
pressure  which,  as  in  equilibrium  calculations,  is  only  expanded  to 
first  order).    With  these  substitutions,  it  is  found  that  the 
fluctuation-dissipation  theorem  is  preserved  up  to  and  including 
o(k  ).    The  resulting  equations  are 

ft  za  "  aijki  3iT7  za  +  Hab(k-;a)  zb  +  >VV 

-  <zb(a1)zc(a2)» 

■  Ra(k,t)  (A>10) 

where 

Hab(*;a)  ■  "^OO  +  *2*%M  ♦  *™M 
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13 

0 

0 

F22 

f23 

0 

0 

f31 

r32 

r33 

L^;(k)  =  r..         r  p  ( A.  1 1 ) 


where  AQ  is  the  thermal  conductivity,  vQ  =  nB/pQ  is  the  kinematic 
viscosity,  V(J  -  |  vQ  +  p~1  Kq  and 


m 

The  vertices  are 


Vabc(k-'21.22)  "  5(^ara2)(Aabc(k.ai)  *  Aacb(kfflg)) 
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"26abea"    (^)ka  3=1 '2  b'c=3"5 

W*-a>  "     "  5  •i*)«!i)«ib)(a>vi0)(J£-a)     -=3-5  b,c=3-5 

0  b,c  <  3. 

(A. 13) 

and  where  we  have  introduced  a  change  of  variables  y  -*•  w  ♦  z.  The 
change  y  +  w  is  defined  by  replacing  the  momentum  density,  j  ,  with 
the  local  velocity  fluid,  u  ,  and  the  total  energy  density,  e,  with 
the  internal  energy  density,  e,  where 

A  a  a 

J    =  p  u 
a  a 

*      1    ^  ^2  A 

e  =  -  p  u    +  £.  (A.14) 

We  note  that 

<6>  =  e0  =  2  P01  J0  +  e0  (A-15) 

so, 
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<V  =  U0a  =  2  *  L  +  <J/P  "  J0/P0> 

<e>  -  e£  =  eQ  +  <e-eQ>.  (A. 16) 

The  second  change  of  variables,  w  +  z,  is  defined  by 

z1  "  cl(p  "  P0} 


Z2  =  C2(e  ~  e0} 


:3  +  i  =  e^)(k)(ua  -  uQa)  (A.17) 


where  {e  1  (k ) }  is  a  set  of  three  mutually  orthogonal  unit  vectors 


with  e      (k)  ■  k  -  k/k  and 


2  2  ?  3Pn 

V'V,'    '<V2>  -^la0 


2       -2  3p0  ,         °0  r  "t 

'"P°  ^ '<o  "  pj  ^  '  Vp  ^ 


2  °°  3eo'po'w;-  u-,8) 


In  these  equations,  „T  -  p"1  (3p0/3T0)/p  Is  the  thermal  expansion 
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coefficient,  cp  is  the  constant  pressure  specific  heat,  cQ  is  the 
sound  velocity  and  hQ  is  the  enthalpy  density.    For  later  reference, 
we  also  define 


°2  9T0  ,  C2  rT-1  h0 


°3  "  °1   3p0    £0  =  C1   Cp0aT  "  p0°v] 


9T0 


where  cy  is  the  constant  volume  specific  heat  and  Y  =  c  /c 

P  u 

Finally,  the  fluctuating  forces  appearing  in  (A. 10)  have 
correlations, 


<R  (k,t)>  =  0 
a  — 


<B2(k,t)R2(k',f)>  -  (2,)3(o*  2  kgT^yk2  «<k*k')«(t-f) 


<».(!t.t)«b(k'.t.»  -  (2rt3CaBT0(2Vl3<1)3  .  UB  -f  nB)«ab] 


x  k2  5(k+k»)5(t-t')  (A. 20) 


where  a,b  =  3-5  and  all  other  cumulants  are  zero  (i.e.,  <<o(k2)). 


APPENDIX  B 
GENERALIZED  EIGENVALUE  PROBLEM 


The  generalized  eigenvalue  problem  is 


We  solve  this  equation  perturbatively  recalling  that 


2 

ck  >>  vk    -  a. 


Thus,  we  expect  the  eigenvectors  and  eigenvalues 


as 


»«"  -  *  »{■»  ♦  ,2  »«•>  .  .... 


The  lowest  order  equation  is 


H     r(oi)  -  ,Wr(a) 
lab^Ob       A1  ^Oa 


The  eigenvalues  and  eigenvectors  of  H1  are  found  to  be 
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with 


-  -ic    =  -  A(2) 


^  -  o  -  tf>  -  ^5> 


,(1)      1  , 

=  7T~  p°°1,  hoc2«  V  °»  0) 

Vd  cQ 
,(2)         1  . 

*  "  ~J~  (pOc1»  h0°2»  -°0'  °»  0) 

Vd  cQ 

(3)  1 

*  =  ^  (h0C2'  "  P0CT  °'  °»  0) 


(*0 

*        =  (0,  0,  0,   1,  0) 


(5) 

*       =  (0,  0,  0,  0,  1).  (b.6) 


These  vectors  are,  in  fact,  orthogonal.  Because  of  the  degeneracy,  we 
find  that 

r(D  ,(1) 


r(2)  (2) 


1 


where  A      is  to  be  determined  at  next  order. 
The  second  order  equations  are 


H  r(ct)  -  i(otMa)  +  ,(00.(0)  ^  -2  ,  3  (a) 
H1ab51b    "  X2    50a    +  A1    Sa    +  k    aijki  »T  V 


+  (H2ab  +k"2aH3ab)4b)'  ( B.  8 ) 


Multiplying  through  by  <|/6)  and  summing  on  a,  we  find  immediately  that 


X2      =  X2      "2  (F0  +  ak    rn}  ^B-9) 


where 


V0C2 


°0 


—  (poClc3  ♦  h0c2Cj()  +  V(J  =  ^  DT  +  V(»  (B.10) 


is  the  sound  damping  constant  and  D?  is  the  thermal  diffusivity.  For 
a»3  ■  3-5,  one  obtains 


(B.11) 


Now, 
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So  one  solution  of  eq.  (B.12)  for  a=3  is 


A(3)  =  D 
A1  DT 


A3a  =  V  (B. 13) 

Equation  (B.11)  then  becomes 


lx\a\2  -  vQk2  -  ar22  ♦  akx  JL]A     =  0 

y 

EX1     k    -  vQk    +  akx  —  ]Aq5  -  ar32Aai<  -  0  (B.14) 
which  is  solved  by 


AHU  =  1  \5  -  A(k) 

Ac),  =  0  A      =  1 


k  k 

A(k)  =  Z~1T  tan~1(k  /ki) 

kxkl        y  1 

kl"kx+ky'  (B.15) 


This  is  not  the  general  solution  to  (B.14)  and  is  clearly  not  defined 
at  kx  =  0.    However,  it  is  possible  to  show  that  the  Green's  function 
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is  independent  of  the  particular  solution  to  (B.I  4)  chosen  and  is  well 


defined  even  at  kx  =  0. 
Thus,  to  order  k  , 


A    =  -ick  +  -  r  k2  +  1  a  k  k  /k2 
I  2    0         2        x  y 


»2  - » 


»3  -  DTK2 


2  2 
X,.  -  v„k    -  a  k  k  /k 


l4  V 


x  y 


A5  =  v0k2  (B. 16) 


and,  to  lowest  order, 


M  .(a) 

5       -  *  a  =  1-3 


(4) 

5        =  (0,  0,  0,  1,  A(k)) 


(5) 

5       -  (0,  0,  0,  0,  1)  (B>17) 


and  the  set  of  vectors  bi-orthonormal  to  {^{a) } ,  {na}  where 
a    ^a      =    a8'  13       lowest  order, 


(a)  ,(a) 
=  9 


(4) 


=  (0,  0,  0,  1,  0) 


(5) 


(0,  0,  0,  -A,  1). 


a  =  1-3 


(B.18) 


APPENDIX  C 
DYNAMIC  STRUCTURE  FACTOR 


The  dynamic  structure  factor,  S(k,u) ,  was  defined  in  eq.  (4.47) 

to  be 


dk  dk, 
~1  —i 

(2tt)6 
dk.dk. 


s(k'^  =  /!.  dt  /  — ]— |  Z(k-k,)Z(k+k2)C pp(k1,t;k2,0)elwt 


2Re  11  dt  /  S(k-k  )5(k+k ,)C    (k  1ft;k?,0)eia)t, 

(2ir)  * 


(C.1) 


Using  eqs.  (4.36)  and  (4.39),  this  may  be  written  in  the  form 


.  dk. 

S(k.«)  =  2Re  /  — Ls^-k  )c  a)  (a)(k> 

(2-ir)-5  iii  -i 

=  a^   S  3)(!i'^  (C.2) 

(a) 

with  o       the  solution  to 


Huu(a,(va'-vu^(a)^^  = 


1j 


<2  5(JcH£1)i^a>(k1)oa1(^l«) 


(C.3) 
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(a) 

Multiplying  (C.3)  by  ^  (k)o(k-k1 ) ,  integrating  over  kj  and 
integrating  by  parts  the  gradient  term  yields 


dk 

J  o2(k-k  ){[-ia)  ♦  A(a)(k  ;a)  -  \  a.  k.  ^-]a(a)  (k,  ,co) 

(2ir)3  1  2    ij  i  3k  -1 


where  we  have  defined         0(a)  (k      ;u>)  =  oCk-j^  )  a(a)(k  ;u)  and  we  have 
used  the  fact  that  the  shear  modes  do  not  contribute  so  the 
eigenvectors  appearing  here  are  _k- independent.    Equation  (C.4)  can  be 
written  as 


3 

/  ^3  o  +  A(a)(k+ki;a)  -  ^  a.J(k+k1)i  ^-]a(a)  (k^  ,„) 


n^c^^k+k^a)}  =  0.  (c.5) 


In  terms  of  o(a)(k,u>),  S(a)(k,a))  is 


S(a)(k,c)  =  J-S-^j,  j  -(a)(k 

(  2tt)  1  --1 


Now,  the  wavelength  of  the  light  is  usually  much  less  than  the  length 
of  the  illuminated  volume  so,  in  (C.5)  and  (C.6), 


k  »  k 


1  (C.7) 
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yielding  the  approximate  equations, 


c(a),  .  -(a),,  . 
S      (k,a))  =  o  (k;u)) 


(C.8) 


The  right  side  of  (C.8)  is  the  same  for  a  -  1  and  2  so 


S(1)(k,u>)  =  S(2)(k,-a,)  =  SB(k,u>), 


APPENDIX  D 
ELIMINATION  OF  THE  THREE  POINT  FUNCTION 


Upon  making  the  factorization  approximation  in  eq.  (6.22),  it  is 
clear  that  only  the  6=1  term  contributes.    Returning  to  eq.  (6.21)  and 
using  the  adjoint  of  the  Liouville  operator  we  find  that 

<(£+  ya(k))  6yg(k))>ss  =  <ya(k)(£_  6yg(k))>ss  (D.1) 

and  for  6=1 ,  we  have 

«  ^  „ 

£_  fiy^k)  =  (£J06yi(k)  ♦  UJlnt6y^k)  +  (£.)ext«y1(k)  (d.2) 


where 


(£_)06yi(k)  =  -  ikj5y2+J(k) 
(£_)6y1 (k)  =  0 

^Jext^!^  "  °'  (D.3) 


So, 
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3 

<(£+  yo(k))«y1(k)>  =  -i    z    k   g        (k)  (d.U) 

j=1    J  J 


and,  from  eq.  (6.12) , 


La1(k,k';a)  =  La1(k;a)  6(k-k')  (D.5) 


where 


3 

La1(k;a)  =  -i    Z    kj  ga  2+j(k)/gl1(k).  (D.6) 


So,  finally,  using  the  explicit  form  of  g    (k),  eq.  (6.27), 


Ln  (k;a)  =  L    (k;a)  =  0 


La1(ii;a)  =  "  nib"  {p0kBT0  ka-2  "  mn0Uat3*y  +  (D'7) 

We  also  note  here  that  the  source  in  eq.  (6.30)  is  easy  to 
evaluate  to  lowest  order  in  the  shear  if  we  note  that 


£ext  "  0  cx-1 


=  2A(D[y2(k)  +  (|  kBT)yi(k)]  a=2 


=  A(r)  a=3  (D.8) 
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and  because  A(r)<xa,  we  can  set  a=0  in  all  other  quantities,  including 
the  distribution  function,  when  calculating  the  lowest  order 
contribution  in  the  shear  rate.  Thus, 


S1e  =  ° 


S2g  =  2<A(r)y2(k)fiy3(-k)>eq  ♦  f^XAtr^  (k)«ye(!<)>eq  ♦  o(a2) 

Sa6  *  <A(r)6yo(^fiyB(-k.)>eq  +  ^ '  °=3~5  (D.9) 

where,  as  indicated,  the  averages  are  now  done  in  equilibrium  and 
yo0O  is  evaluated  at  a=0  for  all  a.  Now,  to  lowest  order  in  the 
shear, 


so  the  angular  averages  require  that 


S28  =° 


S*B  "  "  a<Zj  PjxPjy/Zj  Pj2>eq(6a36^  +  4ai»a«3>  (D'11> 


to  first  order  in  the  shear.    This  is  easily  evaluated  to  give 


Sa6  =  -  (4  m  N  KBT)  a  (6^  ♦  )  af6=3_5  (D>12) 
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and,  the  contribution  to  the  hydrodynamio  matrix,  is 

aB  ' 

-  -  (1  m2»  kBT)  a  (So3Sol)  *  Vl4<,3)(g~\e  *  0<a2)-  (I>-,3) 

But,  to  lowest  order  in  the  shear,  (g~1 )      =  5      for  a  =  3  or  H. 
Thus, 

H  fi    -  o(a  ).  (D. 1 4) 


BIBLIOGRAPHY 


1.  L.  Onsager,  Phys.  Rev.  37,  405  (1931  );  L.  Onsager,  Phys.  Rev.  38, 
2265  (193D;  L.  Onsager  and  S.  Machlup,  Phys.  Rev.  91,  1505 
(1953);  S.  Machlup  and  L.  Onsager,  Phys.  Rev.  91_,  151  2  (1953). 

2.  J. P.  Hansen  and  I.R.  McDonald,  Theory  of  Simple  Liquids  (Academic 
Press,  New  York,  1976);  J. P.  Boon  and  S.  Yip,  Molecular 
Hydrodynamics  (McGraw-Hill,  New  York,  1980). 

3.  R.D.  Mountain,  Rev.  Mod.  Phys.  38,  205  (1966). 

4.  A.-M.  Tremblay,  M.  Arai  and  E.  Siggia,  Phys.  Rev.  A  23 ,  1451 
(1981);  A.-M.  Tremblay,  in  Recent  Developments  in  Nonequilibrium 
Thermodynamics,  Vol.  199  of  Lecture  Notes  in  Physics.  edlteTby" 
J.  Casas-Vasquez,  D.  Jou  and  G.  Lebon  (Springer,  Berlin,  1984). 

5.  D.N.  Zubarev  and  V.G.  Morozov,  Physica  120 A,  411  (1983);  V.G 
Morozov,  Physica  126A,  443  (1984). 

6.  A.  Lees  and  S.  Edwards,  J.  Phys.  C  5 ,  1921  (1972);  E.  Gosling,  I. 
McDonald  and  K.  Singer,  Mol.  Phys.  26,  1475  (1973);  W.  Ashurst 
and  W.  Hoover,  Theor.  Chem.        1  (1975);  T.  Naitoh  and  S.  Ono, 

7nyS;.if  J:Mx'  2,48  (1978);  T*  Nait°h  and  S-  Ono'  J- Chem-  phys. 

'ut  4bit>  (1979). 

7.  D.  Evans,  Mol.  Phys.  37,  1745  (1979);  D.  Evans,  Phys.  Lett.  74A, 
229  (1979);  D.  Evans,  Physica  (Utrecht)  103A,  343  (1980)-  D  

Sqwiq^*  ^  A~'  290  (1980);  D'  J'  Chern'  P^.'  78, 

3297  (19B2);  D.  Evans,  Physica  (Utrecht)  1 18A,  51  (1983)-  B  — 

Holian  and  D.  Evans,  J.  Chem.  Phys.  78,  5147  (1983);  D.  Evans  and 

G.  Morriss,  Phys.  Rev.  Lett.  5±,  1776  (1983) ;  D.  Evans  and  G. 

Morriss,  Comput.  Phys.  Rep.  J_,  299  (1984). 

8.  W.  Hoover,  D.  Evans,  R.  Hickman,  A.  Ladd,  W.  Ashurst  and  B. 
Moran,  Phys.  Rev.  A  22,  1690  (1980);  W.  Hoover,  A.  Ladd  and  B. 
Moran,  Phys.  Rev.  Lett.  48,  1 81 8  (1982);  W.  Hoover,  Physica 

^chU  118A,  m  (1983);  W.  Hoover,  Annu.  Rev.  Phys!  Chem.  34, 


9, 


3D987e(19S80)"  Kim'  C'  M°ntr0Se  and  T'  Litovitz,  J.  Chem.  Phys.  73, 


146 


1  47 


10.  D.  Evans  and  H.  Hanley,  Phys.  Rev.  A  20 ,  1648  (1979);  D.  Evans 
and  H.  Hanley,  Physica  (Utrecht)  108A,  567  (1981);  D.  Evans  and 
H.  Hanley,  Phys.  Lett.  79A,  178  (1980);  D.  Evans,  Phys.  Lett. 
101  A,  100  (1984). 

11.  C.  Trozzi  and  G.  Clccotti,  Phys.  Rev.  A  29 ,  916  (1984). 

12.  J.  Erpenbeck,  Physica  (Utrecht)  118A,  144  (1983);  J.  Erpenbeck, 
Phys.  Rev.  Lett.  52A,  1333  (1984). 

13.  S.  Hess,  Journal  De  Physique  C3  46,  191  (1985);  S.  Hess, 
International  Journal  of  Physics  6_,  657  (1985). 

14.  A.S.  Lodge,  Elastic  Liquids  (Academic  Press,  London,  1964);  J.H. 
Simmons,  R.K.  Mohr  and  C.J.  Montrose,  J.  Appl.  Phvs.  53  4075 
(1982).  '  — 

15.  B.J.  Ackerson  and  N.A.  Clark,  Physica  11 8A,  221  (1983)  and 
references  therein. 

16.  K.  Kawasaki  and  J.  Gunton,  Phys.  Rev.  A  8,  2048  (1973)-  T.  Yamada 
and  K.  Kawasaki,  Prog.  Theor.  Phys.  53,  111  (1975);  D.  Oxtoby,  J. 
Chem.  Phys.  62,  1463  (1975);  M.  Ernst?  B.  Cichocki,  J.  Sharma  and 
H.  van  Beijeren,  J.  Stat.  Phys.  j_8,  237  (1978). 

17.  A.  Onuki  and  K.  Kawasaki,  Ann.  Phys.  (N.Y.)  121 ,  456  (1979);  A. 
Onuki  and  K.  Kawasaki,  Prog.  Theor.  Phys.  Suppl.  64,  436  (1978); 
A.  Onuki  and  K.  Kawasaki,  Physica  (Utrecht)  1 11A,  607  (1982);  A.' 
Onuki  and  K.  Kawasaki,  Ann.  Phys.  (N.Y.)  131~  217  (1981);  A.' 
Onuki  and  K.  Kawasaki,  Prog.  Theor.  Phys.  63,  122  (1980);'  A.' 
Onuki  and  K.  Kawasaki,  Prog.  Theor.  Phys.  69,  146  (1980)';  A.' 
Onuki  and  K.  Kawasaki,  Prog.  Theor.  Phys.  71_,  16  (1984).' 

R.  Zwanzig,  Proc.  Natl.  Acad.  Sci.  U.S.A.  78,.  3296  (1981  ). 

J.  Machta,  I.  Oppenheim  and  I.  Procaccia,  Phys.  Rev.  Lett.  42, 

I  op  (lllV;±"aC*ta'  l'  °PPe^eim  and  I.  Procaccia,  PhysT" Rev. 

II  If,'  361  (1982)  ^        lm  °Ppenheim'  Physica  (Utrecht) 

20'    42  ^fnSftn?'  ?'G;?-,C°hen  and  J'R-  Dorfman'  Ph^.  Rev.  Lett. 
42,  862  (1980);  T    Kirkpatrick,  E.G.D.  Cohen  and  J.R.  Dorfman, 
Phys.  Rev.  Lett    44,  472  (1980);  T.  Kirkpatrick,  E.G.D.  Cohen  and 
J.R.  Dorfman,  Phys.  Rev.  A  22,  950  (1982);  J.  Dorfman  and  T. 
Kirkpatrick,  in  Systems  Far  from  Equilibrium,  edited  by  L 
Garrido  (Springer,  New  York,  1980). 

21.     L.  Garcia-Colin  and  R.  Velasco,  Phys.  Rev.  A  26,  2187  (1982)-  D 
Sahoo  and  A.  Sood,  Phys.  Rev.  A  30,  2802  (19847. 


18. 
19. 


148 


22.    M.  Marchetti  and  J.  Dufty,  J.  Stat.  Phys.  32,  255  (1983);  R. 

Rodriguez,  E.  Salinas-Rodriguez  and  J.  Dufty,  J.  Stat.  Phys.  32, 
279  (1983);  J.  Dufty,  Phys.  Rev.  Lett.  51,  2159  (1983);  J.  Dufty, 
Phys.  Rev.  A  30 ,  622  (1984);  J.  Dufty,  Phys.  Rev.  A  30,  1465 
(1984).  — 


23. 


S.  Hess,  Phys.  Rev.  A.  22,  2844  (1980);  I.  Pardowitz  and  S.  Hess, 
Physica  (Utrecht)  100A,  540  (1980);  S.  Hess  and  H.  Hanley,  Phys. 
Rev.  A  25 ,  1801  (1982);  S.  Hess,  Physica  (Utrecht)  118A,  79 
(1983).   " 


24 


T.  Keyes,  Phys.  Rev.  A  23.,  277  (1981  ). 

25.  D.  Ronis,  Phys.  Rev.  A  29 ,  1453  (1984). 

26.  H.  van  Beijeren,  Phys.  Lett.  105A,  191  (1984). 

27.  D.J.  Evans  and  G.P.  Morriss,  Phys.  Rev.  Lett.  26,  2172  (1986). 

28.  P.  Drazin  and  W.  Reid,  Hydrodynamic  Stability  (Cambridge 
University,  London,  1984);  k.M.  Case,  Phys.  Fluids  3,  143  (1960) 

29.  A.  Einstein,  Ann.  d.  Phys ik  17,  549  (1905);  A.  Einstein,  Annal. 
d.  Physik  19,  371  (1906). 


30. 


L.  Landau  and  E.  Lifshitz,  in  Statistical  Physics.  Vol.  1 
(Pergamon  Press,  Oxford,  1 980)~  ~   


31  *    ^7ynaz^rl?r°S'  D*  Bed6aUX  and  P*  Mazur'  Physica  (Utrecht)  110A, 
1 ^ '  U982);  V.  Morozov,  Physica  (Utrecht)  126A,  443  (1984).   

32.    M.  Van  Kampen,  Stochastic  Processes  in  Physics  and  Chemistry 
(North-Holland,  Amsterdam,  1981). 

33'  iVmq7«;  Math; Phys-    1993  (1978);  R- Fox'  Phys-  Rep-  c  48, 

wy  (1978);  R.  Fox,  Physica  (Utrecht)  112A,  505  (1982).  — 

34.    I.  de  Schepper  and  E.G.D.  Cohen,  J.  Stat.  Phys.  27,  223  (1982). 

35'    Handbook  of  Chemistry  and  Physloa,  edited  by  R.C.  Wenst  (CRC 
Press,  Boca  Rotan,  Fla.,  1982).  v 

36.  A. A.  van  Well  and  L.A.  de  Graff,  Phys.  Rev.  A  32 ,  2396  (1985). 

37.  See,  for  example,  L.  Letamendia,  J. P.  Chabrat,  G.  Nonchi  R 
Kouch  and  C.  Vancamps,  Phys.  Rev.  A  24,  1574  (1981). 

38.  See,  for  example,  D.  Beysens,  M.  Gbadamassi  and  B.  Moncef-Bonanz 
Phys.  Rev.  A  28 ,  2491   (1983).  sonanz, 


1  49 

39.  L.E.  Reichl,  A  Modern  Course  in  Statistical  Physics  (Univ.  of 
Texas,  Austin,  1980). 

40.  J.W.  Dufty  and  R.F.  Rodriguez,  J.  Stat.  Phys.  33,  261  (1983). 

41.  H.  Grabert,  P.  Hanggi  and  P.  Talkner,  J.  Stat.  Phys.  22,  537 
(1980).  — 

42.  R.  Zwanzig,  Phys.  Rev.  124,  983  (1961). 

43.  J.H.  Ferziger  and  H.G.  Kaper,  Mathematical  Theory  of  Transport 
Processes  in  Gases  (North-Holland.  Amsterdam,  1972). 

44.  L.  Landau  and  G.  Placzek,  Physik  Z.  Sowjet  Union  5,  172  (1934). 

45.  J.W.  Dufty,  J.J.  Brey  and  M.C.  Marchetti,  Phys.  Rev.  A  33.  4307 
(1986).  ^ 

46.  D.  Henderson  and  E.W.  Grundke,  J.  Chem.  Phys.  63,  601  (1975). 

47.  T.R.  Kirkpatrick  and  J.C.  Nieuwoudt,  Phys.  Rev.  Letts.  56,  885 


BIOGRAPHICAL  SKETCH 

James  Francis  Lutsko  was  born  on  March  9,  196-4,  in  Tampa, 
Florida.    There  he  received  his  primary  and  secondary  education.  In 
1978,  he  entered  the  University  of  South  Florida.    In  1980,  he 
transferred  to  the  University  of  Florida  from  which  he  received 
bachelor's  degrees  in  physics  and  mathematics  in  the  spring  of  1982. 
In  the  fall  of  that  year,  he  began  graduate  studies  in  physics  at  the 
University  of  Florida. 


150 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it 
conforms  to  acceptable  standards  of  scholarly  presentation  and  is 
fully  adequate,  in  scope  and  quality,  as  a  dissertation  for  the  degree 
of  Doctor  of  Philosophy. 


J.Wl  Dufty,  Chairman 
Professor  of  Physics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it 
conforms  to  acceptable  standards  of  scholarly  presentation  and  is 
fully  adequate,  in  scope  and  quality,  as  a  dissertation  for  the  degree 
of  Doctor  of  Philosophy. 


P .  Kumar 
Associate  Professor  of  Physics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it 
conforms  to  acceptable  standards  of  scholarly  presentation  and  is 
fully  adequate,  in  scope  and  quality,  as  a  dissertation  for  the  degree 
of  Doctor  of  Philosophy. 


N.S.  Sullivan 
Professor  of  Physics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it 
conforms  to  acceptable  standards  of  scholarly  presentation  and  is 
fully  adequate,  in  scope  and  quality,  as  a  dissertation  for  the  deer 
of  Doctor  of  Philosophy. 


P,.<Ramond 

Professor  of  Physics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it 
conforms  to  acceptable  standards  of  scholarly  presentation  and  is 
fully  adequate,  in  scope  and  quality,  as  a  dissertation  for  the  degree 
of  Doctor  of  Philosophy. 


Professor  of  Engineering  Sciences 


This  dissertation  was  submitted  to  the  Graduate  Faculty  of  the 
Department  of  Physics  in  the  College  of  Liberal  Arts  and  Sciences  and 
to  the  Graduate  School  and  was  accepted  as  partial  fulfillment  of  the 
requirements  for  the  degree  of  Doctor  of  Philosophy. 


December  1986 


Dean,  Graduate  School 


